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Preface to Release 1.1 of the PAN AIR Theory Document 

The major changes that take place in this release of the PAN AIR Theory 
Document will be found in the technical appendices. These modifications 
reflect the following technology changes and enhancements that have been 
implemented since the initial publication of this document: 

0 Two new design splines, “source design 2" and "doublet forward 
weighted" have been implemented (cf. appendices D and I). 

0 A discontinuous source spline has been implemented for source 

analysis networks with the effect that total source strength on a 
network is more accurately conserved. Accurate conservation of total 
source strength causes the program to solve more reliably problems 
having fairly coarsely panelled curved surfaces (cf. appendix I). 

0 The control point recession vectors have been modified. The most 

important modification has been to ensure that network corner control 
points are recessed along the bisectrix of the two adjacent network 
edges (cf. appendix G). 

0 The algorithm that performs the abutment search has been 

substantially redesigned to improve its efficiency (appendix F). 

0 The assignment of doublet matching conditions at abutment 

intersections has been completely reworked. The new procedure 
implements a graph theoretic interpretation of an abutment 
intersection and is substantially more reliable than the old 
procedure (cf. appendix F). 

0 The computation of leading edge forces has been substantially 

improved by the implementation of correction factors that correct 
some discrepancies arising from the nonuniform convergence of the 
numerical solution as panel density is increased. 

0 Added mass coefficients have been implemented. 

The authors would like to acknowledge the efforts of the following 
individuals who provided significant contributions to the revisions of release 
1 . 1 : 


D. T. Chiang, for his discussions of the automatic abutment search and the 
two new design splines, 

R. T. Medan, for his discussion of edge forces, and 

D. J. Purdon, for his discussions of control point locations and added 
mass coefficients. 

In addition, we would like to thank Kathleen Crites and Bonnie Smith for their 
efforts on the typing and figure preparation of the revisions. 
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4.0 An Overview of PAN AIR 


4.1 Historical Development of Panel Methods 

In this section, we will discuss the features which distinguish PAN AIR 
from earlier, less complex, panel methods. These features are (a) "continuous 
geometry," (b) linear source and quadratic doublet variation, and (c) 
continuity of doublet strength. We will explain how these features make 
PAN AIR more accurate and reliable than previous methods, and discuss briefly 
the manner in which these items are implemented in PAN AIR. 

Virtually every panel method approximates the configuration geometry with 
panels whose planform is a quadrilateral. Thus, if the panels themselves are 
planar, only a small class of configurations (such as cylinders and flat 
wings) can be described without gaps being left between panels. These gaps 
tend to be very small, except for highly twisted surfaces. In subsonic flow, 
the gaps cause little if any numerical error, but in supersonic flow the 
cumulative effect of the gaps is serious, not because of "leakage" of flow 
through the gaps, but because the doublet strength jumps abruptly from a 
non-zero value to zero at a panel edge which does not exactly meet the 
adjacent edge. In PAN AIR, gaps are closed by means of "piecewise flat" 
panels, that is, panels which are comprised of several planar regions. 

Some panel methods use "curved" panels, generally paraboloidal in shape. 
These approximate the configuration surface far more accurately in regions of 
high curvature such as the leading edge of a wing, but necessarily have gaps, 
even though small ones. Thus they are excellent for the analysis, of subsonic 
flow, but not for supersonic flow. 
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As we stated earlier, PAN AIR employs a linear source variation and a 
quadratic doublet variation. That is, the basis Function b. corresponding 
to a source parameter is locally linear, while the basis function 
corresponding to a doublet parameter is locally quadratic. This contrasts 
with earlier, simpler programs in which the doublet and source variations were 
locally constant. 


The reasons behind the "higher order" singularity distributions in PAN AIR 
are discussed in detail in Appendix B.4. Briefly, they are as follows. 
Consider a control point on a panel, and assume the source and doublet 
distributions in the immediate neighborhood of the control point are 
polynomials. Then we show in Appendix B.4 that a source distribution locally 
of the form 


« 2N 

o(C,n ) = E E 

N=1 i=0 
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N 
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2N-i 


(4.1.1) 


^i r^2N+l-i (4.1.2) 

does not induce any perturbation velocity locally. That is, even terms in the 
polynomial source distribution and odd terms in the doublet distribution do 
not generate a local perturbation velocity. So, since we have concluded that 
constant source and doublet strengths are insufficient, the next reasonable 
higher order approximation to use is linear source strength and quadratic 
doublet strength. 


or a doublet distribution 

« 2N+1 n 

u(?,n) • E E a. N 

N=1 i=0 ^ ^ 


4.1-2 






I^HI 

^■1 

IHI 


o 

B 

B 











»— 1 







1 

1 

1 

1 

1 

z 

o 

p 

i° 

QC O' 
LU 

1 

1 


o 
1— 1 
z- 
o 
to 
a: 

LU 

o 
1— 1 
z 
o 
to 

QC 

LU 

1 ^ 1 






S 




3 

o. 

z 

Q_ 

z 

1 s 1 






CO > 

Ijj • 

coSj 



o 

-J 

to 

to 


1 

1 

1 

1 

1 

LU 
O- 
—1 X 
C LU 

a ^ 

_ *— • 

< LU 

cu to 

*“• 

1 

1 

Lu 

<-> 

z 

o 

to 

QC 

Q 

1 . 
CJ 

z 

<5 

Q 

z 

c 

CJ 

z 

o 











LU 

to 

to 


■ 

■ 

■ 

■ 

■ 

1 1 1 LU 

z > 
2 . 


■ 

■ 

a. 

JS 

CO 

z 

to 

CQ 

z 

to 


■ 

■ 


■1 

■1 

z « 





■■ 

LO 

z 

o 

h — 

t ® 

to 

o 

z 

3 >- 

H 

ii 

>- 
O —1 
z z 
< o 

■ 

to 

CO 

z 

3 

to 

z 

o 

z P 

1-1 c 
z a: 

1 

to 

C3 

Z_,0 

■ 

1 

C_) 

1— I 

oc 

J— 

oo 

OJ 

cc 

1— itOoO 
^cruj 

1 Z|~I 
z;:^Q 
03o 

Z CQ 

0£ Z 

«t o 

z 

< 

-J 

a. 

O Q 

>- —1 
—1 LU 

z z 
c <c 

LU O- 

z 


to to 

CO LU 
Z •— 1 
(—1 Q 
3 O 
CQ 

1 

z 

<t 

z 

c 

O- 

h- Z 

O »-t 
Z LU 

z 

o 

CO 

1 

ZlO< 
«t O- 

—I 
cu 

1 

1 

to 

3 
z o 
O -1 










UH 

z 
»— 1 

BOUNDARY 

CONDITION 

•-I Ll 

< — 1 
o <c 

n:i 

•- o 
o z 

LU 

D. LL. 

to o 

i§ 

a: —1 
o u. 
z 

< 3 

iS 

O LU 

z 

f g 

Oi —1 
O u- 
z 

NORMAL 

FLOW 

NORMAL 

FLOW 

NORMAL FL 
IN LEAST 
SQUARES SEr 

POTENT lA 

NORMAL 

FLOW 

< 3 

13 

O Lu 
z 

ARBITRARY 

<t>, 

ARBITRARY 
0 . V<f> 

■ 













1 













DOUBLET 

TYPE 

NONE 

CONSTANT 

CONSTANT 

LINEAR 

LINEAR 

QUADRATIC 

ftOOTH, CUBI 
QUADRATIC 

CONSTANT 

QUADRATIC 

CONTINUOUS 

QUADRATIC 

CONTINUOUS 

QUADRATIC 

CONTINUOUS 

QUADRATIC 








to 






— 

»- 


h- 

1— 

z 

1— 

z 

o 



oc 

cc: 

a: 

QC 

SOURCE 

TYPE 

<c 

J— 

to 

z 

o 

o 

NONE 

< 

to 

z 

o 

o 

c 

h— 

to 

z 

o 

o 

c 

J— 

to 

o 

CJ 

1 — 

oc 

c- 

<t 

z 

O' 

NONE 

to 

z 

o 

o 

c 

LU 
1— < 

<c 

LU 

z 

< 

LU 

z 
►— « 

LU 

z 

HH 







—1 


_J 

-J 










c 


to o 



to LU 

to LU 







o 





=> to 

Z to 



5 

1— 

c 

_l 

»- 

c 

1— 

c 

o 

— J 

o 

1— 

c 

o o 

Z —1 

z o 

1— 1 CQ 

o 

_J 

o 

1— 

c 

iS5 

orr 

INUO 

CEWI 

'LAT 


U- 

LU 

LU 

Lu 

Lu 

CQ 

c 

Lu 

1“ Di 
Z LU 

CQ 

«c 

Lu 



C3 






oc 

c 

Q. 


O o- 
o >* 
z 

oc 

o. 


O Q. 
CJ 


U- 

B 

CM 

B 

CSJ 

CO 

B 

m 

VO 

n 

<x> 

OY 

■ 


H 


H 



H 

Ei 


B 

19 




1962 



00 

CM 

CO 

CO 


to 

VO 


p 


VO 


VO 

rx. 

Cv 

r«- 

Cn. 


r». 


p 

Ul 

>“ 

o> 

OY 

OY 

OY 

OY 

OY 

OY 

OY 

OY 

OY 

OY 














>- 













z 






LU 



QC 

1— CD 





LiJ 







LU 

q: z 



c 

z u. 

C I-I 


o 

t-H 

0 £. 

3 



o 

z 

z 


o 

CQ 

CO 

z 

LU »— « 
CO ^ 
CQ LU 




z 





QC 

a: 

z 

oc 

z z 



ac 







LU O 



a: <1 

_l = • 


O LU 






Q 

CQ CC 


o 

o. 

< LU 


S 5 

z z 

o 
H-. O 

to to 

O -J 
Z CiJ 

< z 

RUBBERT 
(VORTEX L 

Z CO 
< OJ 

i_ <: 

o 

o: 

< 


z 

c 

to 

1“ 

Qi 

LU o 
3 li. 
to 

• LU 

oc ^ 

LU 

o 

z 

< 

o 

z 

z 

<t 

z 

o 

to 

o 

Z LU 
C Z 

►H 

to _l 
QC 

Q 

1 — o 

LU C^ 

to 1 — 
oc o 

o: 

1—1 

<t 

QC ZC. 
O 1— 

to S 

ffi 

o 

</) 

LU 

CQ 

CO o 

QC Z 

a: 

z 

LU Z 
-U CJ 

LU — 1 
_J 

z 

LU 

3: 

LU 

z 

z 

QC 

o 

LU 

z 

o 

oc 

LU 

f 

o 

z 

3Z < 
UJ Z 

Z Q_ 
LU = 

< 

o. 


4.3-1 (11/30/81) 


Figure 4.1 - Historical Overview of Panel Methods 
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Figure 4.2 - Division of panel into subpanels 
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Figure 4.3 - Panel and far field control point 


5.5 Singularity Splines 


In this section we will discuss without details the construction of spline 
matrices for analysis and wake networks. The technical details of the spline 
construction, and all discussion of splines for design networks, will be 
reserved for Appendix I. In figure 5.13, we illustrate the locations of 
source parameters on a source analysis network, and the locations of doublet 
parameters on a doublet analysis or wake network. 

Source parameters on analysis networks are located at panel centers only. 
Doublet parameters on analysis networks are located at panel centers and in 
addition along network edges as illustrated. The value of a source parameter 
is always the value of source strength at the parameter location, and 
similarly for a doublet parameter. The "extra" doublet parameters occur at 
those points at which an "extra" corner control point was stationed because of 
edge matching considerations (see figure 5.7). Doublet parameters are 
required on network edges (while source parameters are not) because of the 
quadratic variation of the approximation to the doublet strength. A quadratic 
variation causes rapid changes in doublet strength which make extrapolation of 
the doublet values from the interior of the network to the edges ill-advised. 
The source strength approximation is only linear, however. Finally, doublet 
parameters are only located on the upstream edge of a wake network. The 
doublet strength on a wake network is defined to be constant in the streamwise 
direction, and thus doublet parameters are only required on one edge in order 
to define the doublet strength on the entire network. 
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5.5.1 The Matrices and 

The outer spline matrices define the source strength and doublet strength 
at certain points on the panel as linear combinations of source and doublet 
parameters in the neighborhood of the panel. While a single doublet outer 
spline matrix has been found satisfactory for all purposes, it has been found 
that two source outer splines matrices are generally required. One of the 
source outer spline matrices helps define a continuous source distribution 
used in post processing applications, where it is essential for processing 
considerations that source strength be a uniquely defined function on a 
network. The other source outer spline matrix helps define a discontinuous 
source distribution used in AIC matrix construction, where it is important 
that the total source strength be accurately measured by the corresponding 
integral of the splined source distribution. 

To be precise, consider the panel and network in figure 5.14. A source 

s 

outer spline matrix B is a 5 x 9 matrix which gives the value of source 
strength at ^9 terms of the source parameters 

|x^ , i = 1,...,9| located at the nine panel centers marked by a circle. The 
matrix B^ is a 9 x 21 matrix giving the values of doublet strength at Pj...,Pg 

in terms of the doublet parameters i=l 2l| located at the 21 panel 

centers marked by an x. Because p is a continuous locally quadratic 

function whereas a is only a locally linear function, p must be defined at 9 

points on a panel by B while a is only defined by 5 points by B^. The 

values of a at the 5 points are called "panel source parameters," while the 

values of y at the 9 points are called "panel doublet parameters." 
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5.5.2 Definition of SPSPL 


C 

The subpanel spline matrices (one source matrix SPSPL and one doublet 
matrix SPSPL^ for each of the eight triangular regions composing the panel) 
each define the coefficients of the polynomial distribution of singularity 
strength on the triangular region as a linear combination of the singularity 
strengths at the panel points P^ mentioned above. Thus, on each triangular 
region, source and doublet strengths o(C‘, n') and w(^', n') are defined in 
terms of local coordinates (5‘, n') 


o(C, n') = Oq + ogC 

+ an' 
n 

p(5', n') = Pq + 

+ p n' 
n 



where the constants a^, o , o^, 

o’ nn 


1 .2 

L nn 
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spline matrices: 


and 
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5.5.3 Construction of B Matrices for Continuous Singularity Distributions 


A B matrix associated with a continuous singularity distribution is 
constructed one row at a time. Each row defines the singularity strength at 
panel corner, edge midpoint, or panel center in terms of surrounding 
singularity parameters. This identical row vector then becomes part of the B 
matrix of each panel which shares the particular grid. point. This insures 
that the value of the singularity strength is identical as one approaches the 
grid point from the interior of any of the panels sharing it. 

The source strength at a panel corner is obtained from the source 
singularity parameters located at the centers of the four panels sharing that 
corner, as illustrated in figure 5.15. The dependence of a on Xj^,...x^ 
is determined by a bilinear fit procedure described in Appendix I.l. 
Essentially, this procedure determines what "bilinear" function (a bilinear 
function in two variables (^,n) is a quadratic function which reduces to a 
linear function for constant C and n) 


f(E, n) = a + bC + Cn dC n 


(5.5.4) 


is determined by the four values x^-, and then sets to be the "value" 
the function takes at that point. By "value", we mean a row vector 
(apa 2 ,a 3 ,a^), so that 


oi - ^a^ a2 a3 a^j 



(5.5.5) 
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J 



regardless of the values of the x^'s. 


Now, finding the row vector that describes the source strength at a panel 

center is very simple, since a source parameter is located there. To obtain a 

matrix for a panel, we assemble the row vectors corresponding to the 5 

grid points. Each row vector has length 4, but by adding zeros each row 

vector expands to length 9. Thus each row vector has one entry from each of 

the 9 source parameters in the neighborhood of the panel. While only four 

parameters lie in the neighborhood of a particular corner point, (cf. figure 

5.14) nine parameters lie in the neighborhood of at least one of the panel 

S 

corners. Collecting the five row vectors, we have the 5 x 9 matrix B , 
which was first introduced by equation (4.2.8). 


Thus, for the panel in figure 5.14, 

0 ** 0**000 
** 0**0000 
gS ^ 000**0**0 

0 000 ** 0 ** 
0000*0000 


(5.5.6) 


where the columns of B^ are arranged according to the integer labels given 
to the source parameters in figure 5.14. Here, an asterisk denotes some 
generally non-zero entry. 
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The outer doublet spline matrix B*^, introduced in equation (4.2.9), is 
similarly constructed row by row. To obtain the row vector describing y at a 
panel corner, a least squares fit is used. As shown in figure 5.16, y(P) is 
obtained by finding the quadratic function u(e,n) which best goes through the 

12 values at the 12 doublet parameter locations in the neighborhood of P in 
a weighted least squares sense (a quadratic function in two variables 

certainly can not go through 12 values exactly). The computation of the 

weights is discussed in Appendix I. 1.2. 4. The quadratic function thus 

obtained (its 6 coefficients are each row vectors of length 12, since they 

depend on the x?) is evaluated at P to obtain y(P). This weighted least squares 
procedure will be described in detail in Appendix 1.5. 

To obtain a row vector defining y at a panel edge midpoint, we again use a 
weighted least squares fit, though this time we only fit to 8 neighboring 
singularity parameters, as illustrated in figure 5.17. If the grid point lies 
near the network edge, a special treatment (which is described in Appendix 
I.l) is used. 

5.5.4 Construction of the Discontinuous Source Outer Spline Matrix 

The discontinuous source outer spline matrix is constructed by means of a 
two stage process. First, a linear source distribution over the whole panel 
is determined in terms of the panel's neighboring source parameters x^ , i = 
1,...,9 by means of a weighted least squares procedure. Second, this 
distribution is evaluated at the five points Pj^, ? 2 , P3, P^, Pg to give 
the dependency of the five "panel source parameters" upon the neighboring source 
parameters x^. 
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It is the first step of this process that ensures that total source 
strength is accurately measured. This accuracy is achieved by the combinati 
of the linear fit and the fact that the panel's own source parameter is 
heavily weighted in the least squares fitting procedure. 

It is appropriate to observe here that although the discontinuous source 
outer spline is not explicitly constrained to be continuous, it is in fact 
very nearly continuous wherever the configuration is sufficiently finely 
panelled that the angle between adjacent panel normals is less than, say, 10 


5.5.5 Construction of SPSPL 

Next, let us consider the method by which the subpanel spline matrices 
define singularity distributions within a panel. In referring to the panel 
illustrated in figure 5.18, we will write for a(P^) and for 
p(Pi). 


Recall that 02 , a^, and ag are defined in terms of neighboring 
source singularity parameters by the matrix We then define 


^5 - 2 (°i + =7 ^ = I ‘' 4 ^’ ‘^8 = 1 (''4 + 

We have now defined at all vertices of all 8 triangular regions, and we 
now define a linear distribution a(t', n')^, i = 1,...,8 on each triangular 
region by specifying it to be the unique linear distribution to attain the 
appropriate values at the 3 vertices of the triangle. 
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in (5.9.9) over the combined surface S^US 2 , while in fact we only 
integrate the expression over The evaluation of the integral over S 2 » 
the edge force, requires the use of some special extrapolation and correction 
techniques. The basic idea is to evaluate the limit in the expression for 
edge force, (cf. ref. 5.2); 

r— 2 

edge force per unit length = (ir/8)6p [lim (viVTl) ] 

x „-^0 

n 

(here. Bp is so edge normal compressibility factor and x^^ is the distance 
from the edge) by evaluating the expression (w/\^) at panel centers near 
the edge. A correction factor is then applied to the result to account for 
some nonuniform convergence effects arising from the fact that PAN AIR does 
not allow w to behave like cVT^ in the neighborhood of the leading edge. 

For more details of the edge force computation, see appendix 0. 
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dt 


Figure 5.19 - Opposite orientations of adjacent networks 



Figure 5.20 Configuration and image 
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Figure 5,21 Pressure coefficient rules 



Figure 5.22 - Surfaces of integration for leading edge force 
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D.l Networks 

A network is an array (with, say, M rows and N columns) of points in space 
which define a portion of the configuration geometry. In addition, source and 
doublet distributions are defined on the network (that is, the network is a 
"composite" network), with singularity parameter locations and spline methods 
determined by the network's "source type" and "doublet type". 

D.1.1 Network Types 

The possible source types are "analysis", "source design 1," "source 
design 2," and "null", while the doublet types are "analysis", "doublet 
forward weighted," "design", "wake 1", "wake 2", and "null". Source and 
doublet analysis networks are used in conjunction with boundary conditions 
defining impermeability. Design networks are used in conjunction with 
"design" boundary conditions, that is, those which specify tangential 
velocity. Note that a "doublet forward weighted" network is really a doublet 
design network. A network of type "null" is used to denote that the source or 
doublet strength is zero; one could equally well use an analysis network in 
conjunction with the uniform boundary condition 

a = 0 

or y = 0 

To model a wake, as described in section B.2, one would generally use a 
doublet wake network in conjunction with a source null network. The boundary 
conditions, which are only imposed at the wake leading edge, specify the 
matching of doublet strength on that edge to the doublet strength at the 
trailing edge of the adjacent wing network(s). In figures D.l through D.3, we 
illustrate the singularity parameter locations corresponding to each of these 
network types. 

D.l. 2 Wake Networks and the Kutta Condition 

Two types of wake networks are available. In wake 1 networks, the doublet 
strength is variable along the leading edge, and constant in the indicially 
perpendicular direction. In wake 2 networks, the doublet strength is constant 
over the entire network. In the example of figure B.9, the wake extending 
behind the wing would generally be modeled with a wake 1 network, while the 
portion of the wake extending back, from the body would be modeled with a wake 
2 network. 

The two types of wake networks have distinct purposes. The wake 1 network 
is PAN AIR'S approach to satisfying the Kutta condition (see below), while the 
purpose of the wake 2 network if to carry over the doublet strength from the 
wing to the plane of symmetry. 

The Kutta condition, which should hold at the trailing edge, is 

A Cp = 0 (D.1.1) 

where Cp is the pressure coefficient. If the freestream direction is the x 
direction, and the freestream has unit magnitude, then (cf. (C.1.5)) for a 
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(D.1.2) 


thin wing, the linear expresssion for aCp is 

‘Cp = -2 ^ 

Now, the boundary conditions on the wake insure doublet continuity from 
the thin wing to the wake. In addition, it follows from section J.ll that the 
zero normal mass flux boundary conditions along the trailing edge of the wing 
insure the continuity of the x-component of the doublet gradient. 

Now, the wake spline is such that the doublet strength is constant in the 
streamwise direction, that is. 


Since the normal mass flux boundary conditions insure matching of the doublet 
x-derivative, we have, in light of (D.1.2), 


aCp 


trailing edge of wing 


(D.1.4) 


Thus for a thin wing, the use of a wake 1 network results in the 
satisfaction of the Kutta condition, using the linear pressure coefficient 
formula. It is therefore natural to use the wake 1 network to satisfy the 
Kutta condition for a thick wing. This is done in PAN AIR, even in the 
absence of a theoretical justification of its validity. 

Wake 2 networks have a purpose which is not related to the Kutta 
condition. In figure B.9, we show a wake 1 network emanating from the wing 
trailing edge. Now, the body is not a lifting surface, and therefore one 
would not in general expect a panel method to require a wake emanating from 
the body. The wake 2 network is required in PAN AIR, however, because in its 
absence the doublet matching boundry conditions on the wake 1 network would 
drive the doublet strength to zero along its inboard edge. 

Because the doublet strength on the wake is constant, the doublet gradient 
is zero, and thus the surface vorticity, n x vw, is zero. This corresponds to 
the physics of the configuration; that is, the body "sheds" no vorticity. 


D.1.3 Indexing 

We now discuss the indexing system used internally in PAN AIR. The user 
specifies an array [CP(I,J)J of panel corner points, where I, 1 < I < M, is 
called the row index, and J, 1 _< J £ N, is called the column index. ~The upper 
surface is defined by an upward pointing unit normal fi whose direction is the 
vector cross product (direction of increasing column index) x (direction of 
increasing row index). In figure D.4, we illustrate a network with n pointing 
up from the paper. The network edges are labeled in counterclockwise fashion 
as shown, and each panel's corner points are similarly labelled in 
counterclockwise fashion. The point CP(1,1) is called the origin of the 
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F.3 Automatic Abutment Search 


In this section we describe the procedure used to identify "pairwise 
abutments." An edge segment S (the line connecting two adjacent boundary mesh 
points) is said to form a pairwise abutment with a network edge E provided the 
end points s_ and s+ of S satisfy d(t±, E) < e where e i^ some 
user specified tolerance distance. Here, distance from a point s to an edge 
E is defined in the usual way, d(s',E) = m^in d(t,^. 

eeE 


The practical implementation of thp d^finit^on requires that one know how 
to compute the distance from a point s = s_ or s+ to an ^dge E . Let 
edge E consist of edge segments T-j connecting points t-j_]^ and 
ti, i = 1, 2, ... n. Then d(s,E) is given by the formula 

d(f,E) = min dCs,T.) (F.3.1) 

l0£n 

where the distance from a point to a line segment is given by 


d(t,T^) 


- fi-i, - ‘i-i> < ° 

|t-t, I (s- 

|(f - t._j) X (tj - t,_j) I j 


u - ‘i-1 


otherwise 
(F.3. 2) 


Having clearly defined the concept of a pairwise abutment of an edge 
segment with an edge, we now describe what is meant by a pairwise abutment of 
an "edge portion" with an edge. First, by "edge portion" we mean a subset P 
of some network edge consisting of contiguous edge segments, ^k* * * 

The edge portion P then, forms a pairwise abutment with E provided each of the 
edge segments S] does. If the situation illustrated in figure F.4 occurs, 
several pairwise abutments of edge portions P, Q, R with edge E will be 
defined. There are, however, limitations on the permissibility of 
configurations of the form of figure F.4. These limitations are noted in the 
User's Manual, sec. B.3.5. 


Clearly the process of determining all pairwise abutments requires a large 
amount of computation for a configuration with many networks. (The amount of 

work is proportional to where is the number of edge segments in 

the configuration.) In PAN AIR this computational effort is reduced by 
avoiding the computation of the distances d(s±^E) whenever the edge of 
which S is a segment is sufficiently far away from E that a pairwise 
abutment is impossible. 
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F.4 Doublet Matching Along Abutments 

The purpose of doublet matching boundary conditions is to ensure that 
equation (F.1.1) holds at every point along an abutment, even though a 
boundary condition of this form is imposed at only a finite number of points. 

In this section we discuss the enforcement of the doublet matching condition 
at control points along the interior of the abutment while the enforcement of 
doublet matching at the ends of an abutment is treated in section F.5. 

The ability of a finite number of boundary conditions to cause doublet 
matching along the full abutment depends directly on the splining techniques 
used to define the doublet strength along network edges. We discuss this 
subject in section I. 1.2. 5, but we will summarize here the results we derive 
there. 

Given any pair of network edges belonging to an abutment, we call the 
first edge a refinement of the other if, at every point where a panel corner 
is located on the second edge, a panel corner is also located on the first 
edge. According to this definition, each network edge in figure F.5 is a 
refinement of the other, while in figure F.6, edge 1 is a refinement of edge 2. 

We show in section I. 1.2. 5 that if an abutment contains edges Ei,..., E„, 
and some edge Eu is a refinement of each of the other (n-1) edges of the 
abutment, then doublet matching can be forced to take place along the entire 
abutment provided it occurs at the endpoints of the abutments, and at the 
panel edge midpoints on edge E|(. In practice, precise doublet matching will 
not occur because PAN AIR uses a "least squares" rather than a differentiable 
edge spline (see section I. 1.2. 5). The extent to which doublet matching fails 
to occur is very small, and has been found experimentally to be negligible. 

The program takes into consideration the above results when assigning one 
edge of an abutment to be the "matching edge", that is, the edge at whose 
panel edge midpoints doublet matching boundary conditions are imposed. Thus, 
when no special considerations intervene, the edge with the densest paneling 
is assigned to be the matching edge. Assuming that the program user has in 
fact provided one edge in the abutment which is a refinement of all the other 
edges, then that edge is clearly the most densely paneled edge, and so doublet 
matching wi 1 1 occur. 

Under certain circumstances, the program does not assign the most densely 
paneled edge in the abutment as the matching edge. The first such case arises 
from a matching edge of a doublet design or doublet wake network taking part 
in the abutment. 

Unlike doublet analysis networks, design and wake networks are asymmetric; 
boundary conditions are only imposed along certain edges of these networks, 
called matching edges, as illustrated in figures D.2 and D.3. When a matching 
edge of a design or wake network belongs to an abutment, the program assigns 
it to be the matching edge for the abutment, even if it is not the most 
densely paneled edge. This is done mostly for convenience; the user is not 
likely to know what boundary condition to impose at the control points along 
the matching edge and so the program evades this dilemma by assigning doublet 
matching boundary conditions there. 
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The other circumstances under which the most densely paneled edge is not 
chosen as the matching edge is illustrated in figure F.7. To be specific 
whenever the curve defined by an abutment is "supersonic" (that is, no point 
on the edge is in the domain of dependence of any other point), then that 
network edge which is a leading edge (that is, upstream of the remainder of 
the network) is assigned as a matching edge. 

dam I5d is largely empirical. Experience with the 

PAN AIR pilot code with the configuration shown in figure F.8a, illustrates 
the need for imposing doublet matching on the leading supersonic edges of 
networks in supersonic flow. When doublet matching was imposed along the 
trailing edges of networks 1, 2, and 3, the solution was completely erratic 
while shifting the matching boundary conditions to the leading edgL of 
networks 4, 5, and 6 resulted in a solution which was physically reasonable. 

The reasons for the numerical problems resulting from the assignment of 
matching edges as shown in figure F.8a are not precisely known. It is known 
however, that specification of normal mass flux in a two-dimensional 

planar, supersonic flow problem is equivalent to specification of 
the doublet gradient. We may see this by combining equation (C.1.5), which 
states that aCp is proportional to a^/ax, with equations (11-1) and (11-3) 

0 reference F.l, which states that aCp is proportional to normal mass flux. 

Thus for a two-dimensional configuration, the specification of zero normal 
mass flux at panel center points, in combination with doublet matching at the 
trailing edge, is equivalent to the situation in figure C.IO, with the 
trailing edge boundary condition becoming specification of y. But this set of 
boundary conditions, in conjuction with the doublet analysis spline, does not 
have a unique solution. We see this by noting that the doublet distribution 
wo(x) shown in figure F.8b satisfies y = 0 at the leading and trailing 
edges, and ay/ax = 0 at panel centers. Thus, if some solution y(x) exists 
which satisfies the boundary conditions above, so does y(x) + a yo(x) for 
all real numbers a. Thus it is not permissible to specify y at the trailing 

edge of a two-dimensional network on which normal mass flux is specified at 
panel centers. ai. 

Of course, the configuration in figure F.8a is not a two-dimensional one. 
Nevertheless, it seems to have enough resemblance to a two-dimensional 
configuration that the imposition of doublet strength specification on the 
trailing edge of networks 1, 2, and 3 is an unstable boundary condition 
specification for the doublet analysis network spline in use on those networks. 

Summarizing, PAN AIR selects the network edge to be used for doublet 
matching along an abutment according to the following criteria: 


(i) Matching edge of a doublet design network 

(ii) Matching edge of a doublet wake network 

(iii) If neither (i) nor (ii) occur, and the abutment is supersonic the 
leading edge of the most "downstream pointing" network is used 

(iv) If none of the above occur, the most densely panelled network edae 

IS selected for matching. ^ 
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Thus, to insure precise doublet matching, a program user must be sure, for 
every abutment containing a matching edge of a doublet design or wake network, 
that this edge is a refinement of all the other network edges. Similarly, if 
the abutment contains a supersonic edge, the leading edge of the most 
downstream pointing network must be a refinement of the others. Finally, in 
all other cases, some edge must be a refinement of all the others (recall 
that, if two edges have identical paneling, each is a refinement of the other). 

If these rules are followed, the edge chosen by the program as the 
matching edge will in fact always be the most densely paneled one, so precise 
doublet matching will occur. This does not necessarily mean that minor 
violations of the rules will be serious. For instance, in figure F.9 
(ignoring the gap - filling panels for the moment), the edge of the network A 
is not quite a refinement of the edge of network B. There is no reason to 
believe, however, that the doublet discontinuities which result from the small 
discrepancies in figure F.9 are significant. 

Next we must discuss the complications introduced into the above procedure 
by considerations of symmetry. Fortunately, these are few and simple. First, 
we must recognize that either all of an abutment lies on a plane of symmetry 
or else no portion of it lies on a plane of symmetry - an abutment cannot 
partially abut a plane of symmetry. If a network edge lies on a plane of 
symmetry along part of its length and then breaks away, PAN AIR will recognize 
two abutments and place an extra control point at the network course grid 
point at which the breakaway takes place. 

Now when an abutment lies on a plane of symmetry, doublet matching along 
that abutment takes place automatically whenever the potential is symmetric 
with respect to that plane of symmetry. Consequently we find in PAN AIR that 
abutment doublet matching conditions are imposed only on selected symmetry 
conditions when the abutment lies on a plane of symmetry. These are given; 

0 If an abutment lies on the 1st plane of symmetry, impose edge 
matching on only. 

0 If an abutment lies on the 2nd plane of symmetry, impose edge 
matching on only. 

0 If an abutment lies on both planes of symmetry impose edge matching 
on only. 

(Remark: The various symmetrized potentials, etc., are defined 

as follows. The superscripts S or A indicate whether the given function 
is symmetric or antisymmetric in a particular plane of symmetry. The first 
(second) superscript indicates the function's symmetry property with respect 
to the first (second) plane of symmetry.) 
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F.5 Abutment Intersections 


Within the interior of an abutment, the equation 

Esi Pi = 0 (F.5.1) 

can easily be imposed by assigning a particular edge as the matching edge, and 
imposing (F.5.1) at the panel edge midpoints on this edge. 

At abutment intersections, points where two or more abutments meet, (see 
figure F.IO), the choice of points at which to impose (F.5.1) becomes more 
difficult. Only one matching boundary condition may be imposed at a network 
corner point (since only one control point is located there), yet the corner 
point lies at the end of two distinct abutments. We will say that a corner 
point C is "assigned" to an abutment A if the boundary condition imposed at C 
is doublet matching across A. 

A second complication is the danger of overspecification. Consider the 
abutment intersection, formed by four networks, illustrated in figure F.IO. 

Let us define to be the doublet strength at the corner of network 
at this intersection. In order to obtain doublet matching, we require 

Pi = w2 = = 1^4 (F.5. 2) 

But these are only three equations. Thus, if we assign corner point 
to abutment Ai, C2 to A2, and C3 to A3, that is, impose the boundary 
conditions 

ui = U2 ^1 

P2 = v*3 sT ^2 (F.5. 3) 

and 113 = P4 at C3 

we have satisfied (F.5. 2). If we were to assign corner point C4 to abutment 
A4 in addition, the resulting boundary condition 

114 = at C4 (F.5. 4) 

would be redundant, since it follows from (F.5. 3). If a row of the AIC matrix 
corresponding to (F.5. 4) were generated, the resulting matrix would therefore 
be singular, since this row would be a linear combination of three rows 
corresponding to (F.5. 3). 

Thus overspecification must be avoided if the program is to provide a 
numerical solution to the potential flow problem. This is straightforward for 
any reasonable example, but clearly the program must follow a well-defined 
method which assigns corner points to abutments in such a manner that doublet 
matching occurs at all abutments while no overspecification occurs. As an 
example, the abutment intersection in figure F.ll may arise from a realistic 
airplane configuration, yet an automatic procedure assigning corner points to 
abutments is not obvious. In this section, then, we will describe a graph 
theoretic interpretation of this abutment intersection problem together with 
the corresponding solution of this problem. This will be accomplished in two 
phases. In section F.5.1 we will describe the graphical representation of an 
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abutment intersection in the "usual case" together with the corresponding 
abutment assignment procedure. Following this, in section F.5.2 we will 
outline those special features supported by PAN AIR that affect abutment 
assignment together with the modifications to the basic assignment procedure 
that enable PAN AIR to correctly implement those special features. 


F.5.1 Graphical Representation of an Abutment Intersection 

In figures F.IO, F.ll and F.12 we present diagrams for three examples of 
abutment intersections. We will denote by P^j (the abutment intersection 
point), the point at which the various abutments meet. The directed graph G 
associated with the abutment intersection is constructed as follows. 

Let a small sphere S be constructed with P/^j as its center. The nodes 
of G are to be identified with the points at which the various abutments 
pierce S. The branches of G are to be identified with the lines on S along 
which the various networks involved in the abutment intersection cut the 
surface of S. An orientation (direction) for a branch/line is induced in a 
natural way by the orientation of the network that generates it. To see how 
this is done, let N be a network that is involved in the abutment intersection 
and let L|\| (= NHS) be the line along which N cuts S. Denote by N' the 
subsurface of N that lies outside of S. Notice that the line L(\| is part of 
the boundary of N'. Now since N' is a subsurface of N, the orientation of N 
provides an orientation for N'. An orientation for N' in turn provides an 
orientation (that is, a direction of traversal) for the boundary of N'. This 
traversal direction is, of course, the usual counterclockwise traversal of the 
boundary when the network is viewed from above. If then, one traces the 
boundary of N' in the traversal direction provided by its orientation, part of 
the trace will move along the line L|^ in a unique direction. This direction 
is the orientation of L|^. If we denote the abutment at the beginning of the 

line L|^ by Aj^ and the abutment at the end by A,J, then we say that the branch 

induced by network N points from node/abutment A" to node/abutment AjJj. 

The procedure given above generates a directed graph G lying on the sphere 
S. Such a graph can be spread out on a plane simply by puncturing S at some 
point that does not lie on the graph and then stretching the surface S, with 
the graph G imbedded in it, out onto a plane. 

The graphs generated by the abutment intersections of figures F.IO, F.ll 
and F.12 are drawn in figure F.13. 

An alternative (and consistent) interpretation of branch orientation is 
possible if, as is usually the case, the abutment intersection point Pni 
coincides with a corner point of network N lying at the last point of edge 

kj^ and at the first point of edge When this happens, we denote by 


* Remark: Usually it will happen that kjjj = k* where k* = mod (kj^, 4) + 1. 

However if edge k* is a collapsed edge of network N, we will have 
kj, = mod (k,;f + 1,4)+ 1. 


F.5-2 (11/30/81) 



AjiJ the abutment in which edge k|^ participates and by Ajjj the abutment in which 

kiJj participates. Having done this, we say that branch/network N proceeds from 
node/abutment A^ to node/abutment 

Once the directed graph representing an abutment intersection has been 
constructed, it is quite an easy matter to write down all of the doublet 
matching conditions of the form (F.5.1) associated with the abutments in the 
abutment intersection. Given a node/abutment A, we compute the values s^ 
associated with the doublet matching condition according to the following 
rules: 


Si = 


+1 if branch/network Ni is directed away from A 

-1 if branch/network N-j is directed toward A 

0 if branch/network N, is not connected to A 


(F.5.5) 


Using these rules, together with the graphs provided by figure F.13, we obtain 
the following sets of matching conditions 
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(F.5.8) 


(The reader is urged to verify the correctness of these matching conditions by 
carefully re-examining the original figures). 

The matching conditions given above have been written down in a format such 
that they can be readily re-expressed in the shorthand form 
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A(G) w = 0 


(F.5.9) 


where u is the vector of doublet. values associated with the branches/networks of 
the graph and A(G) is called the incidence matrix associated with a directed graph 
G. For the three graphs given in figure F.13, the incidence matrices are 
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Fig. F.13c 
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Having expressed the doublet matching conditions at an abutment intersection in 
terms of an incidence matrix A(G) for the directed graph G describing the abutment 
intersection, we are now in a position to avail ourselves of the many powerful 
results from graph theory. In fact, graph theory not only provides theorems that 
yield much information about the structure of abutment intersections, it also 
provide a number of powerful algorithms that, when suitably tailored, generate the 
doublet-matching boundary condition assignments required by PAN AIR. The standard 
reference for all graph theoretical results quoted in this appendix will be ref. 

Tu ’ u .f*^oory with Applications to Engineering and Computer Science " 

Throughout the remainder of this appendix, it is assumed that the reader has at 
least a nodding familiarity with the elements of graph theory. 

7 2 and 9^6^^ref^^F^2)'^°''' 9>^3ph theory that we shall need is given (cf. THEOREMS 

Theorem* Let G be a connected directed graph containing n nodes and having 

incidence matrix A(G). Then rank(A(G)) = n-1. Furthermore, any set of 
n-1 rows selected from A(G) is a linearly independent set. 


* For the reader familiar with electrical circuit theory, this result is 
equivalent to the result that for a connected circuit with n nodes the 
Kirchoff current law provides n relations, any n-1 of which are 1 inearl v 
independent. 
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This result can be extended to directed graphs G that are not connected by 
observing that any such graph can be written as the union ^connected 
components. Thus, for a graph G with k components, we write 


The theorem can now be applied individually to each component G-j. In fact, 
whenever G is not connected, PAN AIR performs doublet-matching assignments by 
treating separately each component G^j of G. Consequently, to simplify the 
discussion, we shall always assume in what follows that the graph G associated 
with an abutment intersection is a connected graph. 

We now turn to the problem of assigning doublet matching conditions 
(nodes) to replace user specified boundary conditions at control points lying 
on a network (branch) involved in the abutment intersection. At this point we 
treat just the simple "usual case" characterized by the following conditions: 

(i) the graph G associated with the abutment intersection is 
connected, 

(ii) each branch/network in the graph has a doublet distribution, 

(iii) each branch /network in the graph has a control point available 
for use as a doublet matching control point whose hypothetical 
location is essentially coincident with the abutment 
intersection point P/^j. 

When these assumptions are made, the following procedure ensures that n-1 
abutment matching conditions are selected to replace user specified boundary 
conditions on n-1 networks 

A.l Form a spanning tree Tc=G. T must contain all of G's nodes but have 
no loops. (Remark: By a theorem of graph theory, T will contain 

(n-1) branches.) 

A. 2 Select any node of G(T) and label it as the ground node . 

A. 3 Defoliate T by removing one branch at a time until all (n-1) branches 
have been removed. Algebraically, this is accomplished by 
constructing a sequence of trees T = T^ =>T|r,_]^=> ... => Ti where 
Ti is obtained from Tj+i by finding a node of degree 1 in Tj+i 
(not the ground node), and removing that node and the single branch 
to which it is attached. As this is done, the doublet matching 
condition associated with the node is assigned to replace a user 
boundary condition on the network associated with the branch. 

In figure F.14 we illustrate the application of this procedure to the graph 
given by figure F.13b. 
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F.5.2 Modifications to the Abutment Assignment Procedure 

The many special features supported by PAN AIR, especially those features 
associated with symmetry, add considerable complication to the basic algorithm 
for performing doublet matching at abutment intersections. In addition to 
symmetry, the most significant complicating features are, 

(i) an edge of a network may be marked "no doublet edge matching" bv the 
user. 

(ii) a network's corner control point may be a "non-matching corner point" 
in the sense that there is no boundary condition (i.e., AIC row) 
associated with the corner control point that might be replaced with 
a doublet matching condition. 

This presentation of our response to these complications will consist of three 
parts. First, we will describe some of the general considerations that must 
be taken account of. Second, we will outline the numerous ways by which 
symmetry, "no-doublet edge matching" and "non-matching corner points" affect 
the properties of networks, corner control points, user boundary conditions, 
edges, abutments and abutment intersections. Third we will describe the 
algorithm employed by PAN AIR that produces a consistent set of doublet 
matching assignments while satisfying the constraints imposed by PAN AIR's 
special program features. 

General Considerations 

In addressing the problem of symmetry, we adopt the fundamental point of 
view of formulating separately the boundary value problems for the various 
symmetric and antisymmetric parts of the perturbation potential. Thus, when 

two planes of geometric symmetry are present, four symmetrized potentials 

JAA^ JSA^ see appendix K for definitions) will generally be required 
and a separate boundary value problem will be formulated for each. If an abutment 
intersection lies away from any plane of symmetry, doublet matching assignments 
will be the same for all symmetry conditions and indeed will be the same as if no 
symmetry were present at all. On the other hand, if the abutment intersection 
point lies on a plane of symmetry, doublet matching assignments will be performed 
separately for each symmetry condition and it becomes important to know the 
following facts about the abutment intersection: 

0 which plane(s) of symmetry the abutment intersection point lies on, 

0 which plane(s) of symmetry the individual abutments lie on, 

0 which (if any) plane of symmetry an individual network may lie in. 

(N.B. A network is said to lie in a plane of symmetry if all ils points 
lie on the plane of symmetry so tKat the network normal is parallel to 
the plane of symmetry normal.) 

The first of these facts must be known in order to determine which planes of 
symmetry are active , in the following sense: the first (second) plane of symmetry 

IS said to be active if P^j lies on the first (second) plane of symmetry and the 

symmetry condition under consideration is either or $^A second 

plane of symmetry: or P^). This information is important for two reasons. 
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First for an abutment lying in an active plane of symmetry, doublet matching is 
automkically satisfied by virtue of the symmetry properties of the symmetrized 
potentials. Consequently, a doublet matching condition is never explicitly imposed 
for an abutment lying on an active plane of symmetry. Second, the doublet 
distribution on a network lying in an active plane of symmetry is identically zero 
so that such networks do not participate at all in the assignment of doublet 
matching conditions. 

Other Significant Considerations 

Here we outline the ways in which the properties of networks, corner control 
points, etc. are modified by PAN AIR program features. In addition we will discuss 
briefly the mechanisms by which modifications take place. 

Properties of Networks 

0 A network may lie in a plane of symmetry. A network will lie in a plane 
of symmetry if (i) the user explicitly informs PAN AIR of this fact in the 
program input, or (ii) by examination of the network's mesh points PAN AIR 
determines that all points on the network are closer than the geometric 
tolerance distance to a plane of symmetry. 

0 The program user will assign to each network one of the following network 
doublet types: 

Doublet Analysis (DA) 

Doublet Design 1 (DDl) 

Doublet Forward Weighted (DFW) 

Doublet Wake 1 (DWl) 

Doublet Wake 2 (DW2) 

No Doublet (NOD) 

Only those networks not marked "no doublet" are of interest when one is 
analyzing doublet matching at an abutment intersection. 

0 If a network's doublet type is not "no doublet", the user may specify the 
doublet value at all control points. Even if this specification is y = 0, 
such networks are treated differently from "no doublet" networks in 
the sense that these networks are always involved in doublet matching 
along abutments and at abutment intersections. In particular, the 
abutment intersection processing may replace a specified doublet 
boundary condition with a doublet matching condition. 

Properties of Corner Control Points 

0 A corner control point may be a "matching corner point" in the sense 
that no user specified boundary condition is available to generate an 
AIC row that has been reserved for the control point. Thus, the 
control point must have a doublet matching condition assigned to it. 

A "matching corner point" is any corner control point of the 
following types: 

(i) any corner control point lying on the matching edge of 

a DWl network. This includes corners 1 and 2 as well 
as any extra control points along the matching edge. 
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the matching corner control point of a DW2 network 


(ii) 

(iii) any corner control point lying on edges 1 and 4, the 
"doublet matching" edges, of a DDl network. 

(iv) an extra control point on a DA, DDl or DFW network as 
long as that edge is not marked "no doublet edge 
matching." 

0 A corner control point may be a "no matching corner point" in the 
sense that no AIC row has been reserved for any boundary condition 
associated with that control point. Control points of this type 
include: 

(i) corner points 3 and 4 of DWl networks 

(ii) corner points 2, 3 and 4 of DW2 networks 

(iii) corner point 3 of DDl networks 

(iv) any extra control point on an edge that is marked with 
either "no doublet edge matching" (a user 
specification) or "non-matching edge" (edges 2, 3 and 
4 of DWl networks; all edges of DW2 networks and edges 
2 and 3 of DDl networks 

(v) a regular corner control point (that is, one that has 
an AIC row reserved for it) for which both adjacent 
edges are marked "no doublet edge matching." 

0 A corner control point is said to lie on a plane of symmetry if its 
hypothetical location lies on a plane of symmetry. It is possible 
for a control point to lie on two planes of symmetry. 

0 A corner control point is said to lie jii a plane of synmetry if it 
lies on a plane of symmetry and, in addition, the panel normal at the 
control point is parallel to the plane of symmetry's normal. (In PAN 
AIR, the only control points that lie in a plane of symmetry are the 
control points on networks that themselves lie in a plane of 
symmetry. ) 

0 A corner control point always has associated with it two abutment 
ends . (Note: the initial end of an abutment is denoted Fy (+1 ) 

(Abutment index) while the terminal end of an abutment is denoted by 
(-1) (Abutment index).) The orientation of the network provides an 
ordering for these abutment ends, e.g. <A^, A^>, where Aj and 

^2 are abutment end indices. In figure F.15 we illustrate the 

anomalous situation in which A^ = A 2 . In the process of constructing 

the graph associated with an abutment intersection, such a situation 
gives rise to a self-loop , which is subsequently ignored. 

0 The single corner control point at the end of a smooth abutment has 
associated with it two non-smooth abutments as illustrated in figure 
F.16. 
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Properties of User Boundary Conditions 

0 When a doublet matching condition is assigned to a network lying in 
an inactive plane of symmetry, the doublet matching condition will 
replace the user specified antisymmetric boundary condition. An 
antisymmetric boundary condition has the general form 

• ^)a eg u tp . Vu = b (F.5.14) 

(Note: A doublet network lying in a plane of symmetry must be 

assigned an antisymmetric boundary condition of the form given above.) 

Properties of Edges 

0 A network edge may be marked "no doublet edge matching" by the 
program user. When this is done, the network's doublet strength 
along that edge will not participate in any doublet matching 
conditions for the abutment(s) in which that edge is involved. 

0 A network edge may be marked "closure edge" by the user. If in 

addition the user has specified that the closure boundary condition 
override doublet matching, the network's control points will not be 
used for doublet matching along the interior of the closure edge. It 
is generally permissible for control points at the ends of a closure 
edge to be used for doublet matching. 

0 Certain network edges are implicitly marked "matching edge" by PAN 
AIR. These include: 

(i) Edge 1 of a DWl network 

(ii) Edges 1 and 4 of a DDl network 

0 Certain network edges are implicitly marked "non-matching edge" by 
PAN AIR. These include: 

(i) Edges 2, 3 and 4 of a DWl network 

(ii) All edges of a DW2 network 

(iii) Edges 2 and 3 of a DDl network 

0 Every portion of a network edge is involved in exactly one abutment. 
Properties of Abutments 

0 Abutments involving any interior edge of a network are forbidden 

0 The initial end of an abutment may participate in an abutment 

intersection with the terminal end of an abutment. (Figure F.17 
illustrates how this situation can arise for a tube panelled as one 
network.) Since doublet matching must, in general be imposed at both 
ends of an abutment it is necessary to distinguish the initial and 
final ends of an abutment during abutment intersection analysis. As 
noted above, the scheme used by PAN AIR labels the initial end with 
(+1) (abutment index) and the terminal end with (-1) (abutment index). 
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0 An abutment may lie on 0, 1 or 2 planes of symmetry. If any portion 
of an abutment lies on a plane of symmetry, the whole of the abutment 
lies on that plane of symmetry. If an abutment lies on an active 
plane of symmetry, doublet matching along that abutment is 
automatically satisifed by virtue of the symmetry properties of the 
symmetrized potentials Thus, doublet matching conditions are 

never enforced for an abutment lying on an active plane of symmetry. 

0 Smooth abutments do not explicitly enter into the analysis of an 
abutment intersection. 


Properties of Abutment Intersections 


0 An abutment intersection may lie on 0, 1 or 2 planes of symmetry. 
Selection of Matching Conditions 

We now describe the process by which doublet matching assignments are made 
while carefully taking into account the considerations outlined above. 

The basic procedure for performing doublet matching will remain 
essentially the same as the tree defoliation procedure outlined at the end of 
section F.5.1. The complete procedure, however, will be significantly more 
complex at each stage of processing. A summary of the stages of the complete 
procedure is given - 

B.l Construct the graph G containing just those branches corresponding to 
networks whose type is not "no doublet" and do not lie in an active 
plane of symmetry. During this construction process, some 
relabelling of nodes may be performed to account for the "no doublet 
edge matching" feature. 


B.2 


Form a spanning tree for G, Tc=G. 

k 

will have k components, T = U T.- 

i=l ’ 


k 

If G has k components (G = U G^), T 

i=l ^ 

(*). In contrast with the earlier 


algorithm, T is not an arbitrary spanning tree of 6, but rather is 
constructed with careful regard for the properties of the 
neworks/branches of G. 


B.3 Without first selecting a ground node, each component tree T^ is 

defoliated, the removal of each branch providing an association of a 
node with a branch. At the end of this process, there will be left 
one node that is not associated with any branch, and this node 
becomes the ground node for the tree. 


* Note, however, that if some component G-j of G contains only one node, 
all of its branches being self loops, then Tj = In practice, self 
loops cause no special difficulty for the doublet matching problem simply 
because the doublet strength on a network that generates a self loop 
"matches itself." 
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B.4 Finally, the assignments are examined to determine if the requisite 
doublet matching conditions have in fact been associated with 
eligible branches. In addition, it is verified that all of those 
branches that must receive a doublet matching condition have in fact 
done so. 

Several remarks about this algorithm are appropriate before proceeding with 
its detailed exposition. First, steps B.l and B.4 are essentially 
deterministic given the specification of the problem. Steps B.2 and B.3, on 
the other hand are substantially heuristic, step B.2 being ambiguous with 
regard to the choice of spanning tree and step B.3 being ambiguous with regard 
to the defoliation strategy. Although the detailed exposition of the 
algorithm will resolve most of these ambiguities, the solutions presented 
should not be regarded as unique. They are, nevertheless, very good. 

B.l The detailed description of an abutment intersection includes the 
following information 

0 Symmetry condition, or 0^^. 

0 A description of the plane(s) of symmetry that the abutment 
intersection point lies on. 

For each network/control point/branch N involved in the abutment intersection, 
the following information is given. 

0 The abutment ends in which N participates at P/\i are given in 

positive sequence <^\, A 2 >. The two edge segments of N that are 
involved in the abutments Ai, A 2 are denoted Ei, E 2 
respectively. Because of the requirement that Ai and Ao be given 
in positive sequence, an oriented traversal of the boundary of N 
would encounter edge segment E^ followed immediately by edge 
segment E 2 . 

0 "No doublet edge matching" information is given for each edge segment 

Ei- 

0 Let Ti denote unit vectors draj^n al^ng edge segments^E-j, pointing 
away from Pat. Let's = + t 2 )/|ti + T 2 I and let Cq be 

the compressibility axis. Then the "compressibility axis downstream 
parameter" a is defined by 

a =5.^0 (F.5.15) 

This parameter is given for each branch. 

0 Each branch is classified as follows; 


0, no doublet network 

2, doublet network, but no corner control point is available 
to enforce doublet matching 
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3, doublet network with regular corner control point, but at 
least one of the edges E-j is marked "no doublet edge 
matching" 

4, doublet network with regular corner control point. 

5, doublet network with "matching" corner control point 

0 The plane of symmetry jiji which N lies, if any, is given. 

For each abutment-end/node involved in the abutment intersection, the 
following information is given 

0 The global abutment-end index 


0 a flag indicating which planes of symmetry the abutment end lies on 

Given this information, the graph G is constructed in two stages. First, 
a graph Gq is constructed using all those branches in classes 2, 3, 4 and 5 
that do not lie in an active plane of symmetry. Next, each branch of Go is 
examined for any "no doublet edge matching" marks on either end. If such a 
mark is found, the corresponding end of the branch is detached from the node 
to which it is attached and a new node is created. The resulting graph is the 
graph G that we seek. Finally all of the newly created nodes together with 
those nodes (abutment-ends) lying on an active plane of symmetry are marked as 
preferred ground nodes." Doublet matching conditions corresponding to 
preferred ground nodes" are never enforced. 


Without loss of generality, we may suppose that the graph generated in 
step B.l is connected. If in fact G is not connected, we simply perform steps 
B.2, B.3 and B.4 separately on each of G's components. 

Any given connected graph will, in general, possess many spanning trees 
and the problem of determining all of the spanning trees of a particular graph 
IS very difficult (cf. ref. F.2, p. 280). The problem of determining a 
spanning tree is fairly simple, however, and the standard algorithm is given 
in (ref. F.2, pp 277-279). Basically this algorithm proceeds by considering 
in some order, each branch of the graph as a candidate for membership in the 
spanning tree. If a branch causes a closed loop, it is rejected, and if it 
does not cause a closed loop it is accepted. 


It is clear that the spanning tree resulting from this procedure depends 
crucially upon the order in which branches are examined. In PAN AIR, branches 
are considered in order of decreasing priority in accordance with the 
following priority scheme: 


P = 


class 2 branches 
class 3 branches 
class 4 branches 
class 5 branches 


(F.5.16) 


Thus, branches with the highest values of p are considered first for potential 
memebership in the spanning tree. Here, a is the "compressibility axis 
downstream parameter" defined above by equation (F.5.15). 
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B.3 In this step of the algorithm, we must select a defoliation scheme for 
the soanninq tree T. If T has n nodes, then there are precisely n distinct 
defoliation schemes that assign nodes to branches, one corresponding to each 
choice of ground node in subalgorithm A. 3 of section F.5.1. Thus, one 
possibility that suggests itself is to examine each of the n possible 
defoliation schemes and determine which one provides the most suitable set of 
assignments. It turns out that such a complex procedure is not necessary and 
that a fairly straightforward modification of subalgorithm A. 3 generates a 
suitable choice of ground node without an exhaustive search(*). 

The modified defoliation algorithm proceeds as follows. As before, we set 

Tn = T and define a sequence of trees Tp = Tp_i => ... =» Jl Lh 

Ti is obtained from Tj+i by identifying in Tj+i a node of degree 1 and 
removing that node and the single branch to which it is attached. If there is 
mor^lhln oL node of degree 1 available, the choice s made by means of the 
following prioritization scheme based upon node type ("preferred ground node 
or not) and branch class (2,3,4 or 5); 

Branch Class 


Node Type 

2 

3 or 4 5 

"preferred ground node" 

5 

3 1 

not a 

“preferred ground node" 

2 

4 D 

Priority classes for 

removing 

a node/branch combination 


Given this prioritization scheme, the node/branch combination in the 
highest priority class is selected for defoliation. In the event that there 
is more than one node/branch combination, the one with the lowest value of p 
(see above) is selected if the node is a "preferred ground node while the 
highest value of p (cf. equation F.5.16) is selected if the node is not a 
"preferred ground node." 


B.4 The algorithm described above will have achieved a successful set of 
doublet matching assignments provided the following conditions are satisfie 


all "matching corner point" branches (branch class 5) are 
included in the spanning tree T. 


a "preferred ground node" is never assigned to a class 5 branch 
(i.e., priority class 1 is never selected.) 

the ground node actually selected for a tree T is a "preferred 
ground node," if any appear in T. 


(*) Note that a completely exhaustive examination of all popible "latching 

assignments possible for a connected graph G would involve (n . s) separate 
cases, where n is the number of nodes in G and s is the number of spanning 
trees. 
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F.6 Gap-Filling Panels 


Whenever a gap whose size is greater than the user-specified tolerance 
distance occurs^between portions of two or more network edges which form an 
abutment, gap-filling panels are defined. The placement of gap-filling panels 
is illustrated for an abutment containing two network edges in figure F.9. 

We now briefly outline the process by which gap-filling panels are 
constructed, first considering the case of two network edges. First, each 
edge in the abutment is "parametrized." That is, each panel 
the edge is assigned a real number t between 0 and 1 inclusive, where ^ is t 
ratio of two distances. The first distance is the sum of the lengths of the 
oanel edges between the starting point of the abutment and the panel corner 
point in^question, while the second distance is the sum of the lengths of all 
the panel edges on a network edge. Thus t represents the 
entire edge length which one has traveled in preceding from the start of the 
abutment to the panel corner point in question. In figure F.9, each pan 
corner point is given a value t. 


Now, 

edges are 
0 = t. 


suppose the distinct values of t which occur for the two network 


I-O’ 


ti. 


where 


< t]^ < . . .< 


tn - 1 


(F.6.1) 


In the example of figure F.9, n = 9, since there are 11 panel corner points 
other than the initial points, and the values t = .80 and t = 1.0 occur twice. 

Now up to n gap-filling panels may be constructed to fill the gap in the 
abutment For each integer i. 1 < i < n, the quadrilateral region with corner 
points lying on the two network edges, with respective parameter values t 
t- 1 and t = ti is examined. If ti is not the parameter value of a 
corner point (Jor instance, .30 is not the parameter value of any cornerpoint 
nn network A in figure F.9), linear interpolation between corner points is 
used to find the point on the panel edge with that 

or more of the four edges of this quadrilateral region have length greater 
than the user-specified tolerance distance e, a gap-filling panel is defined. 

Thus, for the abutment in figure F.6, no gap-filling panels would be 
defined since the gap size is uniformly smaller than e, and so all potential 
gap-??liinnanels have two edges of length less than e. On the other hand, 
if in figure F.9, e were approximately .05 times the abutment length, seven 
aaLfilling panels would be defined, while two potential ones would be 
3??cardeJ Lcause two edges of the panels would be too short. Te reason that 
very small gap-filling panels are never defined is that numerical difficulties 
could occur in measuring the influence of these panels on control POi^s. It 
is not known at this time under which circumstances the resulting doub 
discontinuity might be significant. 

Next, we address the question of how to define the doublet strength on the 
gap-filling panels so that doublet continuity is attained. The edge splines 
Constructed in appendix I assure that, if one edge is a refinement of the 
other, then the doublet strength matches at points on the two network edges 
with the same parameter value t. Thus we want the doublet . 

constant on the panel along the direction perpendicular to the direction of 

the abutment. 


F.6-1 (11/30/81) 


1 We illustrate this in figure F.18^ The four corner points of the panel 

are and Pi+i on edge A and P^ and Pi+i on edge B. We then 
define Mj as the midpoint of the segment P, P^+j, and ' similarly. 

Then, since P-j, Pi+j, and M-j each has the same parameter value as its 
primed counterpart, it also has the same doublet strength, namely p, p-j+i, 
or pi' respectively. Thus, we have defined p at six of the nine panel 
defining points, and we define p at the remaining three panel defining points 
in the natural way. 

This defines p uniquely on the whole panel (see section 5.5). Further, it 
msu^-es continuity of p on the whole edge P, P^+j and the whole edge 
P Pi+ 1 > since these two gap-filling panel edges are subsets of 
ordinary panel edges. Thus the doublet strength on the network edges is 
defined by a single quadratic function in one variable, and therefore agrees 
with the doublet strength on the gap-filling panel edge everywhere, since it 
agrees at three points. 

Finally we must consider the case of three or more network edges meeting 
in an abutment, as illustrated in cross-section in figure F.19. There, 
network Eq is the most densely paneled network edge, and so, if the user has 
followed the paneling rules, it is a refinement of edges Ei and Eo. We 
construct gap-filling panels as illustrated in order to fill the gap in the 
abutment. 

It still remains to be decided how to define the doublet strength on the 
gap filling panels. We want the doublet strength to be continuous across 
edges Ej and E 2 , while at Eg we want 

^"0 - ^^(1) - u(2) = 0 (F.6.2) 

where p ( 2 ) is the doublet strength on the gap-filling panel spanning the gap 
from Eg to E^. On the other hand, the doublet edge splines and matching 
boundary conditions insure 

Wo - tJl - W2 = 0 (F.6.3) 

where p-j is the doublet strength on edge Ej. Thus the specification 

y(l) = PI 

‘'(2) = n (F.6.4) 

insures doublet continuity (in the form (F.1.1)) everywhere. 

Thus the general procedure used by Pan Air for gap-filling panels in 
abutments with three or more network edges is: 

a. choose the most densely paneled edge (which should be a refinement of 
all the other network edges), 

b. define gap-filling panels in the gap between the finest edge Eg 
and each other edge E^, just as if this were an abutment with only 
two edges, and 
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c. define the doublet strength on the gap-filling panels to be equal to 
the doublet strength on the edge E^. 

Note that this procedure works even if there are only two edges in the 
abutment. Further, it insures that the equation 

Es( = 0 (F.6.5) 

is imposed along all network edges. 
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Network 2 


Network 1 



Figure F.l - Impermissible network intersection 
(in cross-section) 



• = region where viscous 
effects dominate 


Figure F.2 - Leading edge vortex (in cross-section) 
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Figure F.12 - Another 


abutment intersection with 4 abutments 
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abutment intersection 


Figure 



F.17 


An abutment intersection in which both the initial 
and terminal ends of an abutment participate in an 
abutment intersection 
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Mi+1 



Edge A 


Figure F. 18 - Gap-filling panel 
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Fiqure F 19 - Gap-filling panels on an abutment with 3 network edges 
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G.O Control Point Locations 


In figure G.l we illustrate the possible location of points on a network, 
called control points, at which boundary conditions are defined. These 
locations are independent of the source type or doublet type network, 

however we will see in Appendix H that meaningful boundary conditions are not 
necessarily imposed at all control points (for instance, on wake networks, 
boundary conditions are only imposed along one edge). 

Note that controls points are illustrated as being located near but not 

ielocu? o^SotentU at a control point that lies directly on a panel edge 
(see sectioned. 11). Later in this section we describe the procedure used to 
recede control points from the panel edge. 

When a network edge is divided into distinct portions belonging to 
seoarate abutments, as illustrated in figure G.2, an extra control POinT, in 
Addition to those in figure G.l, is defined. The same data are computed for 
these control points as for ordinary control points. 

In order to determine its location, a control 
three categories: panel center control points, edge midpoint control ponts, 

SLSel Corner control points. The latter two categories are only defined 
along the network perimeter. The control point is defined by prescribing a 
"hvDothetical location" (the center point, edge midpoint, or corner point at 
which the control point would ideally be located), and a "recession vector" 
which describes the extent to which the control point is receded into 
subpanel from its hypothetical location. 

The size of the "recession vector" has been determined experimentally. 
Rasicanv it has been chosen as small as possible without causing severe 
Sumerical’error wfrefer to figure G.3 in defining the recession vectors. 
There we show an edge control point as Pr and a corner point , 

control points at other points or edge midpoints the procedure is ident 

Panel center control points are only receded very slightly from the center 
point Pg since the doublet distribution is differentiable and the source 
distribution is continuous at Pg ; as a result (see section J. 11) 

Dotential and velocity induced by the singularity distributions 

Koh;^\/pd at Pn It is Still necessary to recede the control point slightly, 

because Influence coefficients can not be computed for a point lying 
TrZZ’y ma sub-panel edge, because the calculations yield singular results 
there. So, we choose the recession vector to be 


(Ps - P9) ^ (P7 - P9) 
^ ” m 


(G.0.1) 
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Edge control points are receded considerably further because of the 
discontmuitip in doublet derivative, surface slope, and source strength 
which occur at network edges. Thus, for the control point located at Pc we 
define its recession vector to be 


R = 6 

^8 - ^5 

<^9 - ^ 5 ) . - 

\%-h 


-k. -k. \ 

(Po - P.)l 



-' 5 I 

Ps-^r 1 

9 ^5| 

- 8 5^1 

where 

r""" <I^8-P5l 

’ 1 ''9 “ ^5 1 ’ 1^1 

-^5|) ff 

|Pl - Pr-I 

(G.0.2) 
> 0 

6 = • 




1 1 5| 


1/10 


otherwise 


(G.0.3) 


(P vector R bisects the angle between the vectors ("Po - Pc) and 

the recession_^vector R emanating from 

point Pg would lie on the line segment joining points Po and Pq. If edge 1 
IS collapsed as in figure G.4, then s is taken to be a tenth. ^The recession 
this°ra^p'*^°^Tf be used for any control points located at Pj and P^o in 
oSssihlfipJ nf .col apsed. then a will be at most a tenth^and 

the panel is skewed. The recession vectors for other edge 

locat ons*^a?rp ^['®.^®tined analogously, control points whose hypothetical 

and s! withdrawn, respectively, into triangles 6, 7, 

edge''?he“ecess"on Sect'or'ir*"* “ collapsed 





This particular construction provides the recession vector ^ with 
properties similar to the edge control point recession vector in (G 0 2) 

Note that the recession vector in (G.0.4) lies in subpanel 1. The recession 
vectors for other corner control points are handled in a similar fashion. 



j 


Some geometric quantities in addition to location and hypothetical 
location are computed by the program for each control point. One of these is 
on which the control point actually lies. This is needed later 
(see section J.8) to insure that an average potential and velocity are 
computed correctly in measuring the influence that the sub-panel on which the 
control point lies exerts on the control point 
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Also, for each edge or corner control point at which a matching boundary 
condition is imposed, a set of "extra hypothetical locations" and their 
associated sign is computed. These arise from the matching boundary condition 
(see section F.l) 

Esi ui =0 (G.0.5) 

where the v-j are the values of doublet strength on different networks. 

In figure G.5, we illustrate an abutment containing three network edges. 
Although the control point is receded from the edge, the matching boundary 
condition involves singularity strengths at the edge; in this example (6.0.5) 
becomes 

2 

Si y(Hi) = 0 (G.0.6) 

i = 0 

where Hq is the (default) hypothetical location of the control point, while 

and H 2 are extra hypothetical locations. 

In Appendix F we indicate how the signs Si are computed; here we 
describe the computation of the extra hypothetical locations. Hypothetical 
locations are computed one abutment at a time by parametrizing the abutment 
(see section F.4), a process that assigns to each panel corner point or edge 
midpoint P on that portion of a network edge belonging to an abutment a real 
number t(P) between 0 and 1. 

In figure 6.6, we illustrate an abutment with two network edges. Given 
the control point and default hypothetical location Hq , we compute the 
extra hypothetical location Hj as follows. Parametrization of the abutment 
gives us t(Ho), and also assigns a value t to every panel corner point and 
edge midpoint on the edge of network 1. 

By interpolation between these points, we find Hj^ as the point satisfying 

t(Hi) = t(Ho) (G-0.7) 

In addition to the coordinates of the extra hypothetical locations, the 
program determines the panel and subpanel on which each extra hypothetical 
location lies, so that the doublet strength can be computed there later. 

Our discussion of matching boundary conditions has assumed we are dealing 
with doublet matching. In the case of a source matching boundary condition 
(which may occur on the edge of a source design network) 

Esi ai = 0 (G.0.8) 

and extra hypothetical locations and signs are computed as before. 
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There are two final pieces of geometric data, associated with control 
points, which we have not yet discussed. These are the normal and conormal of 
the subpanel on which the control point lies. The normal is needed in 
post-processing to compute velocity from the potential and the normal mass 
flux (see Appendix N), while the conormal n is needed to compute the normal 
mass flux from the velocity influence coefficient matrix by the formula 


The computation of the normal is described in section D.2, while 

" " (G.0.10) 

where 

[Bo] - I + (sb2 - 1) Co cj (S.0.11) 

by equation (E.3-.9). 
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I.O 


Singularity Splines 


Singularity splines define the source and doublet distributions on the 
entire configuration in terms of the source and doublet singularity 
parameters. These distributions are defined by a collection of matrices. 

First the source and doublet distributions on a subpanel (recall from section 
5.5 that a panel is partitioned into eight subpanels) are each defined by a 
"subpanel spline matrix" (denoted, respectively, SPSPL^ and pSPL^^) in 
terms of five "panel source parameters" and nine "panel. doublet parameters. 
Thus there are eight of each of these matrices associated with each panel. 

Next, the panel source and doublet parameters are defined by "outer spline 
matrices" (denoted, respectively, by B and B ) in terms of singularity 
parameters located in the neighborhood of the panel. Each panel has 
associated with it a continuous doublet spline matrix, a continuous source 
spline matrix and, possibly, a discontinuous source spline matrix. 

The subpanel spline matrices are defined by equation (5.5.7). That is, 
[SPSPL^]^^^ relates the three coefficients o^, o^, (which define a 
linear source distribution on the subpanel) to the five panel source 
parameters (that is, the source strengths at the five points on the panel 
illustrated in figure I. la). Similarly, [SPSPL^]^^^ relates the six 
coefficients y^,..., y (which define a quadratic doublet distribution on 
the subpanel) to the nine panel doublet parameters (whose locations are 
illustrated in figure I. lb). 


The outer spline matrices are defined by equation (5.5.8). That is, 
rBS] defines the five panel source parameters in terms of the neighboring 
source singularity parameters (generally nine in number), while [B^J defines 
the panel doublet parameters in terms of the neighboring doublet parameters 
(generally 21 in number). 

The subpanel and outer spline matrices are used in the influence 
coefficient calculations. The subpanel spline matrices are first used in 
order to compute "panel influence coefficient" (PIC) matrices (see sections 
4.4.2 and J.l), and the PIC matrices are multiplied by the outer spline 
mkrices to obtain potential and velocity influence coefficient matrices 
([(^IC] and [VIC]) which give the perturbation potential and velocity at a 
point, in terms of all singularity parameters, due to all the panels in the 
configuration (see sections 4.2.3 and 5.9.1). 

In section I.l we discuss the construction of outer spline matrices. 

While their construction is simple in principle, based on a least square 
procedure, in practice it is quite involved because there are many special 
cases. In particular, a special "edge spline" is used near network edges, 
which in conjunction with the doublet matching boundary conditions discussed 
in Appendix F, results in precise matching of doublet strength along network 
edges. In section 1.2 we describe the construction of subpanel spline 
matrices. 
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discuss full panel and half panel spline matrices. 

These matrices define source and doublet distributions specified by sinqle 
the whole panel or half the panel, respLtively, also in 
terms of the panel singularity parameters. The distributions are rouqh 
approximations to the 8 subpanel distributions defined by the subpanel 

calcu^ation^^^ "intermediate field" influence coefficient 

Next, in section 1.4, we discuss "far field moments", matrices describino 
integrals of the singularity strength over a panel in terms of the panel ^ 

singularity parameters. The matrices are used in far field influence 

^he theory of 

We now briefly discuss the reasoning behind the subpanel and outer soline 

stahlp'^^^th" First (cf. section C.3), we require the spline^to be 

disturbance in the singularity distribution caused by a 
perturbation of a boundary condition should die off quickly. 

linearly accurate and the doublet 

spline quadratically accurate. That is, if the source parameters are definpd 
ch source distribution defined by the spline matrices 

th^dlblersol Ls """'“sous property should hold for 

xne aouoiet sp ines. The justification for using a linear source and 

quadratic doublet distribution is given in section B.5. 

Hictlih'"?- spline must be local in nature. That is, the singularity 
distribution on a panel must depend on a reasonably small number of 

0 herl^: problems „h?ch would occur 

?e5u™ed for eaJh’p^nel? would otherwise be 

doublet strength should be continuous (see section 8.4). It 
would be preferable to have continuously differentiable doublet strenqth 

woul^ri^nprm-r^'^f® Strength and smooth geometry as well, since these conditions 
would permit a further integration by parts of the influence coefficient 

Uhtortunately, these go!ls aJfnot 

ach evable without an unacceptable increase in the cost of evaluation of 
influence coefficient integrals. Moreover, while it is a fairly 

to achieve a continuous source distribution it has 
been found that without smooth geometry, continuous source spl ines’induce 

total source strength on a network, seriously 
degrading the accuracy of the aerodynamic influence coefficient matrix"!^ 

Finally, the entries of the PIC matrices, which are defined as sums of 
integrals, should be computable in closed form. That is, numerical 
integration should not be required for the evaluation of the inteqrals The 
reason for th,s reqnirement is one of simplicity. The integrlnd^Jn 

T^‘°- ui k tar too singular to be integrated numerically as thev stand 
It might be possible to partition the integral into a regular part, ^tegr^ble 
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It is the avoidance of numerical integration, in combination with the 
maintenance of geometric continuity, that causes much of the complexity of the 
spline construction. Geometric continuity between panels can be maintained 
either by breaking up a panel into flat subpanels, or by defining a single 
curved panel. The integrals over the curved panel are not computable in 
closed form in supersonic flow, however. 

Once one has decided to use flat subpanels, a minimum of five planar 
regions (those of figure 5.2) is mandated to achieve geometric continuity 
while avoiding any kink in the surface near a panel center control point. The 
use of eight subpanels has been chosen because it offers a convenient method 
of defining a continuous doublet distribution, while not requiring polynomials 
of degree greater than two. An explicit polynomial distribution has been 
chosen rather than a parametric distribution because the integration in 
parametric coordinates can not be performed in closed form. 
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I.l 


Outer Splines 


There are two basic methods that are used for the construction of outer 
splines, corresponding to the cases in which the resulting singularity 
distribution is required to be continuous, or is not. In PAN AIR, a 
continuous outer spline is constructed for the doublet distribution while both 
types of outer spline are constructed for the source distribution. The 
discontinuous source spline is used in the computation of influence 
coefficients and the evaluation of boundary conditions while the continuous 
source spline is used just in the post processing modules. (It was originally 
intended that the continuous source spline be used for all purposes. However 
its inability to conserve total source strength led to the introduction of 
discontinuous source splines in the solution portion of the code* Because the 
post-processing modules had built into them the assumption that source 
strength is single valued, the continuous source splines were retained for 
these essentially less demanding functions.) 

The construction of continuous outer splines is a two step process. In 
the first step, row vectors SP^ and SP^ (called "spline vectors") are 
formulated for grid points as illustrated in Figure I.lc for some typical 
cases. These row vectors define the source or doublet strength at each 
enriched grid point in the network as a linear combination of surrounding 
singularity parameter values. In the second step, matrices [B^J and [B^J 
are constructed for each panel, giving the source or doublet strength at the 
appropriate grid points on the panel (panel singularity values) as linear 
combinations of values of singularity parameters in the neighborhood of the 
panel . 

Thus the matrices [B^] have five rows while the [B^] have nine rows, 
since the source strength is defined at five points on a panel (the panel 
source parameter locations) by row vectors SP^ , while the doublet strength 
is defined at nine points by row vectors SP^ (see Figure I.l). The number 
of columns in a matrix B is variable: it equals the total number of 

distinct singularity parameters on which the panel source or doublet 
parameters depend. The matrices B are assembled from the required row vectors 
SP in a fairly straightforward manner described in the maintenance manual (see 
the preface of SUBROUTINE VECUNM of the DQG module). Briefly, first row 
vectors SP are computed for every grid point in the configuration (except that 
row vectors SP^ are not needed for panel source parameter locations) and 
stored on disk. Then, when the spline quantities for a single panel are being 
computed, the five (or nine) row vectors for each of the panel source 
(or doublet) singularity parameter are fetched from the disk. These row 
vectors are then amalgamated into a single matrix B^ (or B^^) by VECUNM. 

In this section the discussion of continuous outer spline computation will 
simply describe the computation of individual row vectors SP^ or SP^. The 
basic principle is simple. For source splines, the source strength at a grid 
point is fit in a linearly accurate manner to as few surrounding source 
parameters as possible while for doublet splines we do the same in a 
quadratical ly accurate manner. But while the basic principle is simple, 
implementation is complex because of a myriad of special cases which do not 
fit the general rules. 
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The construction of a discontinuous source outer spline matrix is also 
performed in two steps. First, a linear function of the form n + c+ 

^ pJoceduJe that“fSfni5l”il'’uses 

the global source parameters in the neighborhood of the panel 
lron<^ ^ coordinates on the mean panlu) 

linear function is evaluated at the five panel source parameter 
locations, each evaluation generating a row of the source spline matrix, B^. 

1.1.1 Source Splines for Analysis Networks 

1. 1.1.1 Source Spline Vectors for Continuous Splines 

Computing the row vector describing the source strength at the center nf a 
panel in an analysis network is particularly simple, Sinaia LurL 

singularity parameter is located there. Thus the source strength is iust thP 
singularity parameter value; that is, strength is just the 


1 S 

ag = = lIj \ 


or lSP^j = lIj 

the row vector of length 1 with unit value. 


(I. 1.1a) 
(I. 1.1b) 


I I. 1.1. 2 Neighboring Singularity Parameters 

Next, to find the source strength at a panel corner, we perform a 
bilinear fit (a process to be described below) to the four surroundina 
source parameter values. In Figure 1.2, we show the variety of cLer^hich 
may occur in the course of determining the four neighboring source parameter 
locations In the -standard" case (A), the four source pa?aSs a^^the 

choices, and neighboring parameters must be 
obtained by reaching toward the interior of the network. The logi^used fo^ 

points B or C, however, when extended to D, results in a large number of 
nejghbonng source parameters. To keep storage to a minfI„mf„rjto«r(l„ a 
tayly arbitrary manner) from this set of points, those points which are as 
far as possible in index from the uncollapsed edjes of the neLork 

I LI. 1.3 Computation of a Local Coordinate System 

Next we compute the source strength at a panel corner in terms of the four 
surrounding singularity parameter values, once we have in fact located these 
four parameters. The first step is to firm a local coordintu s«tem whSH 
( t, n) plane is the one in which the source strength is to vary linearly. 
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The four singularity parameter locations determine (generally 
non-orthogonal ) basis vectors % and connecting pairs of panel edge 
midpoints, for this coordinate system as follows. Let Pq be the grid point 
at which we wish to find the source strength. Then for any point P on the 
network, we want to be able to determine coordinates 5(P), n(P), t (P), such 
that 


(P-Pq) = C(P) + n(P)% + c(P)v^ ( 1 . 1 . 2 ) 

The (C, n, c) coordinates used here (in section I.l) are not related to the X' 
coordinate system (also denoted ( C, n, c) at times) used in Appendix E aj^d 
secti^on 1.2. Here, ‘v- is a vector perpendicular to Jthe pj^ane spanned by ‘v, 
and “v-. Such a vector is, of cours^, a multiple of "v^ x "v^, but a simple ^ 
dimensional argument shows that if v^ is to be independent of the scale of 
coordinates (that is, if is to be doubled when every point coordinate in 
the network is doubled), we must have 


Vg X Vn 

V, X t, 1/2 


(1.1.3) 


Now, to find the functions which define C, n» and C, let us first take the 
cross product of (1.1.2) with on the lejft, and then the dot product with 

V. Next, we tak^ the cross product with "Vj, on the right, and then take the 
dot product with Since 

X v^ = 0 = % X v„ 

= 0 = v^.v^ 

(v^ xv^ )»v^ = 0 = (v^ x‘v^).'v^ (1.1.4) 

we obtain 

X (P - Pq) 

(P - 'Pq) X 7^ 

(P - Po) • \ 


= n (P) (v^ X V,) 

= C(P) Cv^ X -v^) 

= C(P) (v^.v^) (1.1.5) 


Dotting the first two equations with v^ , we have 
v^ X ("P - Pq) . v^ = n (P) { 7 ^ X IT^) . "v^ 


I 3/2 

n (P) I 7 ^ X \T^| 


( 1 . 1 . 6 ) 


Thus, 


5(P) 


(P - Pq) X V 



n*'' 

nr 


c 


?(P) 



(1.1.7) 
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I. 1.1.4 The Bilinear Fit at Panel Corner Points 

Now, let Pj,...,P^ be the four source parameter locations for any of 
the cases illustrated in figure 1.2, and let = o(P.) denote these 
source singularity parameters. Then these four values of source strength 
define a "bilinear" function in l and n, that is, a function 


'(C, n) = 




( 1 . 1 . 8 ) 


(where the symbol as used in section I.l has a different meaning than in 
section 5 or the remainder of Appendix I), which takes on exactly these four 
values. The function a ( C, n) is defined by the fitting condition 
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( 1 . 1 . 10 ) 

The points Pj,...,P^ are not coplanar in general, and thus computing 
o(C,n) by (1.1.10) in terms of surrounding source parameters, ignores the 
C —component of the parameter locationsj in other words, we project the 
parameter locations to the plane defined by v"^ and "v^. This is justifiable 
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in view of the fact that a reasonable number of panels should be used in 
defining a geometric surface, and thus the distortions due to curvature can be 
neglected locally. 


Now, by (1.1.7), ?(Po) = n(Po) = 0. and thus by (1.1.8), 

o(Pq) = ®o (1.1.11) 

Thus by (1.1.10), 

r S'! 

h 

a(Po) = ao = l 1 0 0 Oj , 

yj 

(1.1.12) 


Setting 


lSPSj = l1 0 0 Oj [bl]-i 

we see that 

fx^l 

a(Po) = ' i ‘ 


(1.1.13) 


(1.1.14) ■ 


Note that lSP^J is just the first row of [BL]"^. Now, by (1.1.14), 

SPS is just the row vector we seek; namely, it gives the value of source 
strength at the point Pq as a linear combination of four neighboring 
singularity parameters. A spline vector may similarly be constructed for 
every panel corner point in the network, whereupon matrices B^ may be 
computed for each panel as discussed at the beginning of section I.l. 

This concludes the discussion of continuous source spline construction for 
source analysis networks. Two special cases, networks with only one row or 
column, and networks with only one panel, are discussed in the maintenance 
document (see section 4-1.4 and SUBROUTINE ONDFIT of the OQG module). 


I. 1.1.5 Discontinuous Source Analysis Splines 


The construction of a panel's discontinuous source analysis spline, 
required for computation of influence coefficients, is achieved by solving for 
coefficients oq, On weighted linear least squares problem 
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where the quadratic form G is defined by 

6 = E t«i(% ^i ^ % ’’i ■ 
i=l 

S 

Here, w^. denotes a weight and is a global source parameter at some point 

in the neighborhood of the panel. The coordinates (Ci, ni) are obtained 
by performing a length preserving projection of ]T-j onto the panel's mean 
plane followed by a transformation of this projected point into the mean plane 
coordinate system. Passing over for the present the selection of w^ and 

x^. , we observe that the minimization problem we have posed has a solution of 


the form 
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Here, we have used the conventional notation A+ to denote a matrix 
pseudo-inverse. Now, letting (Ck, nk), k = 1,...,5 denote mean plane 
coordinates of the projection* of the five panel source parameter locations 
onto the mean plane, we may evaluate the five panel source parameters by 

Ok = ao + C k ^k 

This evaluation process induces an expression for ak in terms of xf, 
providing a definition for [B^], the source outer spline. 



^ ^ = [bS] X 


* This projection is performed in scaled coordinates. 
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We now return to the problem of selection of weights and global source 
parameters to be used in the fitting process. In Figure 1.23 we illustrate 
the variety of cases which may arise when identifying the neighboring source 
parameters to be used in spline construction. For any given panel, the 

parameters are selected as indicated. The weights w-j are chosen as follows: 

c . S • 

if XT does not lie on the panel for which B is 

bein^ computed 

C 

if x^ lies on the panel 

The choice of a very large weight for the panel's own source parameter is 
crucial in that it is this condition that causes total source strength for the 
panel to be correct to sufficient accuracy. 

The foregoing procedure will fail to provide enough data points if the 
network in which the panel lies has only one row or column of panels. When 
this happens, points "pi are selected as indicated in Figure 1.24. The value 

of x^ used for points'?^ that are not global source parameter points is just 
the value of source strength at the panel center. 


1 

10000 


1.1.2 Doublet Spline Vectors for Analysis Networks, Doublet Forward 
Weighted Splines 

Doublet spline vectors SP^ are more complex to compute for a variety of 
reasons. First, the requirement of quadratic accuracy forces the doublet 
strength at a grid point to depend on a greater number of singularity 
parameters than the source strength. Second, to insure doublet continuity 
across network edges on non-smooth abutments (along which boundary conditions 
specifying the matching of doublet strength are imposed), we require that the 
doublet strength at any point on a network edge depend only on the singularity 
parameters located on the network edge. The example of a thin wing with a 
curved planform illustrates the need for this requirement (cf. Figure 1.3). 

The doublet strength is zero at the singularity parameter locations on the 
free network edge. If the doublet strength at a panel corner point on the 
edge depended on singularity parameters in the interior, it could not be zero, 
independent of conditions in the network interior, as we wish it to be. But 
by insisting that it only depend on edge parameters, we insure that it is zero. 

A third cause of increased complexity in determining doublet spline row 
vectors is the introduction of "smooth abutments." These are abutments 
consisting of portions of two distinct network edges, along which splines 
rather than boundary conditions are used to enforce continuity of doublet 
strength. 

For grid points which do not lie on a network edge, obtaining the row 
vector SP^ which describes the doublet strength at each grid point in terms 
of surrounding singularity parameters is a two step process. First, the set 
of surrounding singularity parameters is determined. Second, the doublet 
strength at the grid point is determined in a quadratically accurate manner in 
terms of the neighboring singularity parameters. 
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Doublet forward weighted (DFW) spline vectors are calculated in the same 
manner as the doublet analysis spline vectors. The only difference being 
weighting factors used in the least squares fit. A description of the 
different weighting scheme is given in section I. 1.2.4. 


I. 1.2.1 The Quadratic Least Squares Fit for Panel Corners or Panel Edge Midpoints 

This quadratically accurate procedure is somewhat more complex than the 
bilinear fit employed for source splines. While there is generally a bilinear 
function which exactly fits values at four points (unless three of the four lie on 
a line, which is unlikely if they are panel centers), a quadratic function is less 
well behaved. There is a unique quadratic through six points, unless these points 
all lie on two lines. With very regular paneling, however, it is quite likely 
that six center points chosen as neighbors of a grid point will, in fact, lie on 
two lines. Thus, the procedure we choose for the quadratic fit is a "least 
squares" procedure. 

That is, we choose an excessive number of neighboring singularity parameters, 
and find the quadratic function which takes on the values of the closest 
singularity parameters exactly, while taking on the the values of the remaining 
singularity parameters in a "least squares" sense. The row vector SP^ for the 
grid point is determined by the value the fitting function takes on at the grid 
point, expressed as a linear combination of the neighboring singularity parameter 
values. 


We now describe this least squares procedure more precisely. Let 

(x9, i = l,...,k) be the singularity parameters (in the neighborhood 

of the selected grid point) to which we fit the quadratic function exactly. 

Let (x9, i = k + 1,..., k+m) be the remaining neighboring singularity 
' n 

parameters. Let x^ be located at (C^, c^-)» where the 

computation of these coordinates will be discussed shortly. Once again, 
however, this is not the local (5, n> ^) coordinate system denoted X' in 
Appendix E. 


Let A be the matrix 


i 


"1 . 

1/2 

tl "I 

1/2 

1 

^k 

\ 

1/2 <-1 

^k \ 

1/2 4^ 


This is the matrix for which any function 

f(C, n) = fl + fzi + f3 n + 1/2 f 4 ^2 + fg Cn + 1/2 fe 


(1.1.15) 


(1.1.16) 
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taking on the values 


f( c,. n, 

) = 


i = 

X j • • • 

satisfies 





[^]kx6 

^2 

► 

i 

"2 , 


^6 

V. J 





(1.1.17) 


(1.1.18) 


Now, whenever k < 6 , as it will be in the current applications, equation 
(1.1.18) does not fully specify the coefficients of f. The coefficients are 
completely specified by requiring the minimization of 


I = 


k+m 

E 

i=k+l 


wj [f (C^. n^) - 


(1.1.19) 


where w^ is a "weight", to be discussed shortly, which depends on the 
relative locations of the singularity parameter and the grid point. 

If we write 



"1 

^k+1 nk+1 

1/2 tk+l 

^k+l’ik+1 



[Aj]mx6 = 

_1 

^k+m ^k+m 

1/2 

^k+m’^k+m 

l/2nk.nl. 

(1.1.20) 

equation (1.1.19) 

becomes 





+ 
1— » 

Wi^ 

( ^ i-k, 

V S=1 

\ 

,s)fs “ ^i 

i 



(1.1.21) 


The method by which we minimize (1.1.21), subject to the exact conditions 
(1.1.18), is called a "constrained least squares" procedure, and is discussed 
in section 1.5. The result of performing this procedure is a (6 x (k+m)) 
matrix LSQ such that 
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r 


rx^ 

^1 

• 

= [LSQ]’ 

• 

^6 


k-hn 


( 1 . 1 . 22 ) 

Now, we will construct our (C, n, 5) coordinate system such that at the grid 
point Pq we have 


C(Po) = 0 
n(Po) = 0 


(1.1.23) 


Thus, 


f(Po) - fi 


|^LSQj,.j 


» 

0 

L" k-tm. 


(1.1.24) 


But we required the row vector SPD to define the value at Pq of the 
quadratic function f which satisfied (1.1.18) while minimizing (1.1.21), and 


L$P j = |_LSQj,.j (1.1.25) 

that is, SpD is the first row vector of the matrix LSQ defined by the 
constrained least squares procedure. 

In describing the construction of SP*^ for a grid point in a network 
interior, we have deferred the discussion of three items. These are the 
determination of the set of neighboring singularity parameters, their 
(C, n, ?) coordinates, and the corresponding weights w,-. We will discuss 
them in order as follows. 


I. 1.2. 2 Neighboring Points for Least Squares Fit 

Figure I. 4a illustrates the location of neighboring singularity parameters 
for grid points which do not lie near a network edge. Note that, since a 
singularity parameter is located at each panel center, the spline vector 
SP*^ for a panel center point is (like the spline vector SP^) a vector of 
length 1 with a unit entry. 

Now, if the grid point (panel corner or edge midpoint) lies near (but not 
on) a network edge, the set of neighbors must include singularity parameters 
on the network edge. Actually, we fit a quadratic function to neighboring 
grid points, where these grid point need only be singularity parameter 
locations when they are in the interior of the network. The value of doublet 
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strength at those grid points which lie on the network edge depends in turn on 
a small number of singularity parameters located on the network edge. This is 
a procedure which is defined in detail in the maintenance document (see 
section 4-1). 

Figure I.4.b illustrates the neighboring points we use for the quadratic 
fit to obtain a spline vector for a grid point which lies near (but not on) a 
network edge which does not belong to a smooth abutment. Recall that only one 
singularity parameter is located on a collapsed edge of a network. 

Figure I. 4c illustrates the set of neighboring points when the grid point 
in question lies near a smooth abutment. In this case, we see that the set of 
"neighboring points" may lie in two distinct networks. This is because the 
singularity parameters on the network edges on the smooth abutment (though not 
at the corner points at the ends of the abutment) have been removed for 
reasons of economy. The neighboring points in the same network as the point 
Pq are chosen in the usual manner (see Figure I. 4a) while those in the 
adjacent network are chosen as illustrated. The precise method by which the 
latter points are picked is described in the maintenance document (section 
4-1. 2. 1.2). 


I. 1.2. 3 Construction of a Local Coordinate System 

Next we discuss the construction of a (C, n, C) coordinate system. This 
system is similar but not identical to that of section I. 1.1. 3, and totally 
distinct from the X' coordinate system X'jof Appendix E or section I. 2. 2. 5. 
First, we construct basis vectors w and \T as illustrated in Figure 1.5. 

That is, vV and "v span pairs of enriched grid points adjacent to the base 
point Pq. ^Next, we define t by (1.1.3). Then, analogously to (1.1.17), we 
define 

t(p) = ( p - Pq) " ,(p) = \ " (p - Pq) X(p) = 

|v^xvJ3/2 |7^xv^|3/2 ^ %1 

(1.1.26) 

The bars indicate that these are preliminary coordinate values which will 
be adjusted to account for surface curvature. Consider a cylindrical surface, 
as illustrated in Figure 1.6. If we use the coordinates c and n above, we are 
essentially projecting the surface down to the tangent plane at the point 
Pq. When the surface is highly curved, this makes the points A and D appear 
to be closer to Pq than they really are, since we are dealing with their 
projected images A' and D'. The points B' and C are also closer to the grid 
point than B and C, but not in the same proportion. 

We rectify this by scaling the and n coordinates of a point according to 
its height above the tangent plane. We define a scaling factor 

A(P) = [ C(P)^ ^ n(P)^ ^ C(P)^] 

I(P)2 + n(P)2 (1.1.27) 
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Then we define 

C(P) = A(P) t(P) 

n(P) = A(P) -fr(P) (1.1.28) 

This coordinate scaling assures that the contribution of distant points is 
more accurately measured. Note that the denominator of A(P) is non-zero 
provided P ^ Pq . ' ' 


I. 1.2.4 Weights for the Least Squares Fit 

Next we consider the weights Wi for the least squares procedure. 

In order to provide stability, we would like to fit more closely to nearby 
singularity parameters than further ones. This is done in part by fitting the 
quadratic function exactly to the nearest parameters, as illustrated by Figure 

1 • H • 


A second consideration in determining weights is the desire to give 
heavier weights in supersonic flow to points which are upstream of the grid 
point than to those which are downstream. This weighting has been found 
experimentally to reduce instabilities which arise at high Mach numbers. In 
recognition of these requirements, we set w^ = w(P) where 


w(P) = 1 kMjl - Cq . (P - Pq)/ I P - Pq|) 

2 kMj (1.1.29) 

The constant k is set to zero in subsonic flow in view of the lack of a 
preferred upstream direction. That is, the compressibility direction Cn may 
be replaced by -Cq without changing the solution to the equation. In 
supersonic flow, k has been chosen by experiment, and has order of magnitude 1. 

Since the dot product of unit vectors lies between -1 and 1, the numerator 
of (1.1.29) lies between 1 and (1 + 2k M„ ). Thus the ratio of weight 
(neglecting the effect of distance) at a directly upstream point to that at a 
directly downstream point is 1 + 2k M^ . For Mo<> = 3, and k = 1 (the 
provisional choice for k), this ratio is 7. 

The weights for the doublet forward weighted splines are obtained from 
equation 1.1.29 by setting k = 1 and Mo. = 2. Thus, the simple expedient of 
changing the weights in a least squares fit transforms a doublet analysis 
spline into a doublet design spline (DFW). 


I. 1.2. 5 Edge Splines for Non-Smooth Abutments 

Finally we consider grid points (panel corner points or edge midpoints) 
lying on a network edge. A network edge is divided into distinct portions 
belonging to different abutments. A doublet parameter is located at the grid 
points which form the endpoints of the portion of the edge belonging to the 
abutment (if such an endpoint is not a network corner point, the doublet 
parameter is an "extra" singularity parameter (see figure 5.13)). Doublet 
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pardineters are also located at the panel edge midpoints unless the abutment is 
a smooth one, in which case the parameters are removed (doublet parameters 
located at abutment endpoints are retained for simplicity). 

We discuss first the case of a non-smooth abutment. In that case, the 
value of doublet strength at a grid point depends only on singularity 
parameters located on the network edge. 

Consider the abutment illustrated in Figure I. 7a, with one network edge 
panelled more finely than the other and with a panel corner on the more finely 
paneled network located wherever the more coarsely paneled network edge has a 
panel corner. The goal is to find a splining method such that the imposition 
of doublet matching boundary conditions of some or all of the control points 
on the edge results in the exact matching of doublet strength on the whole 
edge. 

Experimentation with least-squares-type splines shows that they cannot 
satisfy the above considerations. Let us consider, on the other hand, a 
differentiable spline. Let the edge be divided into n intervals, as 
illustrated in Figure I. 7b. It is reasonable to ask how many differentiable 
functions exist, defined by a single quadratic on each of the n intervals. 

Now thoro dre 3n linodrly indGpendent quddrdtic functions . dltogGthGr (sincG d 
quadratic function on an interval has 3 coefficients), and requiring 
continuity at P2.---»Pn-l yields (n-1) constraints on the Set of 
functions, while requiring continuity of derivative at these points provides 
(n-l) additional constraints. Thus, there are (n+2) linearly independent 
piecewise quadratic functions with continuous derivatives. 

But this is equal to the number of control points on the edge, and so 
there is a unique differentiable function which takes on a prescribed set of 
(n + 2) values at the midpoints of the intervals and the endpoints of the edge. 

We can apply this result to the situation illustrated by Figure I. 7a. The 
doublet distribution on the edge 1 will consist of some differentiable 
function defined by a single quadratic on each interval of edge 1. If we now 
impose doublet matching boundary conditions at the control points of edge 2, 
we obtain on edge 2 the unique differentiable doublet distribution defined as 
a single quadratic on each interval of edge 2, which agrees with the doublet 
distribution on edge 1 at the specified points. But, since every interval of 
edge 2 is a subset of a corresponding interval on edge 1, the doublet 
distribution on edge 1 satisfies the above criterion too. So, since the 
distribution is unique, the doublet distributions on edge 1 and edge 2 are 
identical . 

Summarizing, we have shown that if edges 1 and 2 form an abutment, and the 
paneling on edge 2 is a "refinement" of the paneling on edge 1 (that is, every 
corner point of edge 1 is also a corner point on edge 2, though edge 2 may 
have additional corner points), then the imposition of doublet matching at the 
control points of edge 2 results in exact matching of doublet strength along 
the entire abutment. Generalizing, if several network edges meet in an 
abutment, and one edge is a refinement of each of the other edges, then the 
imposition of doublet matching boundary conditions on that edge forces the 
alternating sum of the doublet strengths to zero: 
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where = + 1. 

Unfortunately, the differentiable edge spline, while leading to doublet 
continuity under a greater variety of circumstances, does not permit forward 
weighting in supersonic flow. As a result, the differentiable edge spline is 
insufficiently stable and cannot be used in Pan Air. This fact was determined 
fairly late in the development of Pan Air, and thus a discussion of the 
differentiable edge spline has been included in this document. 

The spline which is actually implemented in Pan Air is a one dimension 
quadratic least squares fit. Consider, for instance, a network edge as 
illustrated in figure I. 7b. The points Pi,..., Pn+U and Mi, 
i=l,...,n, are singularity parameter locations, and the doublet strength 
there is defined to be equal to the value of the singularity located there. 
Thus the doublet spline vector SP^ for each of these points is a unit vector 
of length one, as it is for panel center points in a doublet analysis network. 

Next, the doublet strength at the points P^, i=2,...,n-l, is obtained by 
a constrained least squares analgous to that described in section I. 1.2.1, but 
in one dimension. That is, the quadratic function f(t) (t a variable defining 
distance along the network edge) is found such that 

f( Mk) = w(Mk) = xf k=i,i+l (1.1.31a) 

and f(M^) = y(M^) = r = i-1, i+2 (I. 1.31b) 

in a least squares sense. 

Then the row vector SP'^ which defines w(P-) in terms of the singularity 

D ' 

parameters (Xj^, k= i-1,..., i+2), is such that 


v{P^)=f{P^) (1.1.31c) 

We now discuss the differentiable edge spline which is not implemented in 
Pan Air. First we must compute the spline matrices which correspond to this 
differentiable piecewise quadratic distribution. It can be shown numerically 
that such a function, if it has a non-zero value at one panel center, and is 
zero at all other panel centers and the endpoints of the edge, is never 
identically zero. Rather, it behaves as illustrated in Figure 1.8; 
oscillating with an amplitude which diminishes rapidly but never reaches 
zero. Thus, the spline is stable under doublet specification boundary 
conditions; however, it is not local, since the doublet strength on an 
interval depends weakly on the doublet strength at a panel center far away. 

In order to avoid storing lengthy spline vectors, we must redefine our 
doublet parameters to make the spline local. That is, a doublet parameter on 
a network edge will not have as its value the doublet strength at its 
location. For this purpose, consider the interval [-1, 1] on which we 
define the quadratic function 
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(1.1.32a) 


y(x) = a + bx + cx^ 

Now, y(-l) = a - b + c 

y(l) = a + b + c 

^ (-1) = b - 2c 

dx 

^ (1) = b + 2c 

dx 

Thus, 

n(-l)+^(-l). = a-c = p(l)-^(l) 

dx dx 


(I. 1.32b) 


(1.1.33) 


(1.1.34) 


Generalizing (1.1.34) to the interval [P-j, Pi+i3 in Figure 1. 7b, we 

have 

p(Pi) + 1/2 y(Pi) (Pi+i-Pi) 

p(Pi+l) + 1/2 y(Pi+i) (Pi - Pi+l) (1.1.35) 

We thus define 

x.'^ = y(P.) + 1/2 Vy(P.) . (P.^.^ - P.) 

i = l,...,n (1.1.36) 

x^*^ = y(P^), i = 0, n+1 (1.1.37) 

Now, by (1.1.35) 

X^_l = w(p.) + 1/2 7u(P^) . (Pi_i - P^) 

i = 2,..., n+1 (1.1.38) 

Combining (1.1.40) and (1.1.42), and noting that Vy is continuous, we have 
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(Pi -P,.l) ^(Pi> * 1/2 ^(P,) (Pi*i -P,) (Pi -Pi.i) 


(1.1.39) 


MP,*i -P,) w(P,) -1/2 ^P,) (P, -P,,j) (P,.i-P,) 


Thus, 


Pi - 

Pi-1 

lI 

1 

1 

1— » 

+ 

Pi+i - Pi 


Pi+i - Pi 


J3J .... 1 . 

Pi - Pi-1 + 

hn - Pi 


»? 




(1.1.40) 


This defines the spline vector for Pi as computed with the differentiable 
edge spline, which is not implemented in Pan Air. 


Now, we wish to evaluate v(Mi), i=l,...,n. Again consider a function 
u(x) on [-1, 1] defined by (1.1.31). 

Then 


w(-l) = a - b + c 
u(l) = a + b + c 


w(-l) + ^ (-1) = a - c 

dx 


So, 

p(o) = a = - 

1/4 w(-l) + 1/4 y(l) + 1/2 [y(-D + (-1)] 

dx 

Applying this to Figure I. 7b, we see 


v(M.) = 1/4 w(P.) + 1/4 w(u.^j) + 1/2 x° 

i = l,...,n (1.1.43) 


(1.1.41) 


(1.1.42) 
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This defines the spline vector for M^. 


Equations (1.1.37), (1.1.41), and (1.1.43) together describe u at grid 
points which lie on a network edge belonging to a non-smooth abutment as 
linear combinations of neighboring singularity parameters on the edge. We 
again point out that this procedure is not implemented in Pan Air. 


I. 1.2. 6 Edge Splines for Smooth Abutments 

We now describe the computation of spline vectors for grid points lying on 
smooth abutments. Once again, to obtain matching of doublet strength we 
require that one network be paneled as a refinement of the other, as in the 
example of Figure I. 7a. Then, spline vectors for grid points on the more 
coarsely paneled edge are computed first, followed by spline vectors for grid 
points on the more finely paneled edge. 

Spline vectors for grid points on the coarser edge are also computed by a 
constrained least squares procedure, even though again the "neighboring 
points" lie in two networks. Figure I. 9a shows some representative examples 
which illustrate the procedure for choosing the set of neighboring points. 

The method is described precisely in the Maintenance Document (see Appendix I 
of section 4) . 

Now, continuity of doublet stength along a smooth abutment is insured by 
requiring the doublet strength at a grid point on the more finely paneled 
network to be identical (as a linear combination of surrounding singularity 
parameters) to that at the "corresponding" point on the coarsely paneled 
network. We determine the corresponding point by "parametrizing" the 
abutment, that is, assigning to each grid point a real number t, 0 < t < 1, 
which specifies the proportion of the total abutment length that the grid 
point is distant from the starting point of the abutment. This procedure is 
discussed in more detail in the maintenance document (see SUBROUTINE PRMEDG of 
the DQ6 module). 

Figure I. 9b illustrates the parametrization of an abutment. Now, some 
grid point P'-j on the fine network will have parameter value t'^, where 

tj < t ' -j < tj+]^ ( 1.1.44) 

for some integer j, that is, the corresponding point on the coarse network is 
not a grid point. But, u must vary quadratically on the panel edge, so we can 
obtain p(P'i) as a linear combination of v(Pj). vj(Pj+l)« vi(Pj+2). 

Now, it follows from (1.1.31-32) that on an interval [-1, 1], 
y(x) = a + bx + cx^ = 
y(0) + (l/2y(l) + 1/2 y(-l))x 

+ [1/2 y(l) + 1/2 y(-l) - y(0)]x2 (1.1.45) 
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We can apply this to the interval in Figure I. 9c by making the 
transformation 


^^■"2 - (1.1.46) 
X = 2x' - 1 




t 


j+2 



Equation (1.1.46) maps the interval 
turn is mapped to [-1, 1] by (1.1.47). 


(1.1.47) 

in Figure I. 9c to [0, 1], which in 


Substituting (1.1.47) in (1,1.45), we obtain 
p(t) = y(Pj+i) + [1/2 p(Pj+2) - 1/2 p(Pj)] - tj+2 - t- 


^j+2 


- ti 


+ [1/2 p(Pj+2) + 1/2 u(Pj) - p(Pj+i)] 


2t - tj+2 - t 


'j+2 


“ t -i 


(1.1.48) 


Setting t = t'-j, we have u(Pi) as a linear combination of u(Pi), 
u(Pj+l), and m(Pj+2): ^ 


p(P'i) = (-1/2 +1/2 t2)p(Pj-) 

+ (1 - t 2) y(Pj+^) + (1/2 + 1/2 t 2) y(Pj+2) (1.1.49) 


where T = ^ ^ i ^j+2 ~ 

^^■"2 - tj (1.1.50) 

This concludes our discussion of spline vector construction for doublet 
analysis networks. We have now discussed the computation of doublet spline 
vectors for all enriched grid points in a doublet analysis network. In 
practice, these vectors are all computed and stored on a disk. Then, within a 
loop over panels, the spline vectors corresponding to the nine panel defining 
points are retrieved from the disk, and merged into an outer spline matrix 
by VECUNM. 
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1.1.3 Doublet Spline Vectors for Wake Networks 

Singularity parameter locations for doublet wake networks are illustrated 
in Figure D.3. In addition, if the edge of the wake I networks on which 
singularity parameters are located forms part of more than one abutment, an 
extra singularity parameter is located at the abutment endpoints lying in the 
interior of the edge. 

The purpose of a doublet wake 1 network is to model a wake surface on 
which the doublet strength is constant in the streamwise direction. Thus, 
spline vectors for grid points are constructed as follows. Fii^st, spline 
vectors are constructed for each grid point on the edge containing singularity 
parameters, just as though the edge were part of a non-smooth abutment of an 
analysis network (it should be noted in passing that smooth abutments are only 
permitted between analysis networks). Then, the spline vector constructed for 
a particular grid point on the edge is also used for every grid point lying in 
the column or row of points emanating in an indicially perpendicular direction 
from the edge. This produces a doublet strength which is constant in one 
indicial direction, as desired. In general this direction is the direction of 
increasing row index, though this program default may be overridden by the 
user. See section 7, record N12, of the User's Manual. 

Doublet wake 2 networks are used to define a constant strength doublet 
sheet, whose strength is the value of the one singularity parameter in the 
network. Thus, the identical spline vector is constructed for every grid 
point on the network; namely the row vector of length one with unit entry. 


1.1.4 Source Splines for Design Networks 


I. 1.4.1 Source Design 1 

Only one type of source outer spline, a continuous one, is used for source 
design 1 networks. Singularity locations for source design networks are given 
by Figure D.l. Since a source parameter is located at every panel corner, the 
spline vectors for these grid points are just unit vectors of length 1. 

Spline vectors for panel centers are also straightforward to compute: 

SPS = 1/4 1/4 l/4j (1.1.51) 

That is, the source strength at a panel center is defined as the average 
of the source strengths at all the panel corners. 


I. 1.4. 2 Source Design 2, Discontinuous Source Splines 

For source design 2 networks source parameters are located at those edge 
midpoints on edges parallel to the matching edge. To reduce the complexities 
of splines PAN AIR imposes two restrictions on source design 2 networks: They 
may not have collapsed edges and they may not have just one column or one row 
of panels. 
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The discontinuous source spline for a source design 2 network, used for 
influence coefficient computation and boundary condition evaluation, is 
computed by means of a three stage process. First, spline vectors are 
computed for the panel centers and for those panel edge midpoints that are not 
source parameter locations. (See figure D.lc for an illustration of the 
source parameter location on a source design 2 network.) Second, the five 
source values on the panel, (the panel center and four edge midpoint values) 
are fitted to obtain a source distribution function of the form a + o_ ? + a n. 

0 5 r\ * 

Third, this distribution is evaluated at the five panel source parameter locations. 


This process can be summarized by the equation 


1 



1 0 0 0 0 
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^3 
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(3x5) 
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1 10 0 0 
7 7 


Here, each stage of the process is represented by a matrix. 



The first stage of this process requires further explanation. First the 
source strength at the panel center is taken to be the average of the two 
global source parameters that lie on the boundary of the panel. Next, the 
extra panel edge midpoint source strengths are obtained by means of a bilinear 
fit of neighboring global source parameter data, as illustrated in figure 
1.25. The bilinear fit performed here is essentially the same as the bilinear 
fit described in sections I. 1.1.3 and I. 1.1.4 in connection with continuous 
source analysis splines. 


I. 1.4. 3 Source Design 2, Continuous Source Splines 

The continuous source spline for a source design 2 network, used by PAN 
AIR'S post processing modules for pressure, force and moment calculations, is 
generated by computing spline vectors for the panel center and corners. 

Taking the source strength at the panel center to be the average of the 
panel's two global source parameters, we have 

SPS = l1/2 1/2j 

The panel corner spline vectors are obtained by means of the usual sort of 
bilinear fit using global source parameter data as indicated by figure 1.26. 
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1.1.5 


Doublet Splines for Design Networks 


Figure D.2 shows the location of singularity parameters on a doublet 
design network. For grid points in the interior of the network, spline 
vectors are computed by fitting to neighboring points, as illustrated in 
Figure I. 10a. For grid points on the "matching edges" (which have singularity 
parameters located at the panel edge midpoints), the doublet analysis edge 
spline of section I. 1.2. 5 is used. 

The only unusual aspect of doublet design splines is the edge spline for 
non-matching edges. The doublet parameters are located at panel corners along 
these edges (rather than panel edge midpoints) for stability, since the 
boundary conditions in the vicinity of non-matching edges tend to be doublet 
gradient boundary conditions. For nonmatching edges, as for matching edges, a 
differentiable edge spline and a least squares edge spline are available, 
though once again the least squares spline is implemented in Pan Air. The 
least squares spline is similar to that for matching edges, except that now it 
is at panel edge midpoints that the doublet strength is defined by least 
squaring to the four surrounding edge doublet parameters, while at panel 
corners the doublet strength is defined by a unit spline vector. 

We now discuss the construction of the differentiable edge spline. Let 

be the value of the doublet parameter located as a panel corner as 

illustrated in Figure 1. 10b. We define a row vector ly^j , 0 i 4: n (n the 
number of panel corners on the network edge) of length n as follows. We 
define yq Yn to be row vectors with the entries 1 in the first entry 
and the nth entry, respectively, and otherwise zero. For 1 < i < n - 1, we 
obtain Yi by performing a one-dimensional least squares fit to the 4 (or 3, 
if i=l or n-1) neighboring singularity parameters on the edge. 

Thus, at each edge midpoint, and at the endpoints of the edge, a row 
vector Yi is defined. This is analagous to the situation for the 
differentiable doublet analysis edge spline. We now obtain u at corner points 
and edge midpoints by using the doublet analysis edge spline, but in terms of 
the Yi rather than the singularity parameters. 

For example, we have, analagously to (1.1.40), 


p(P5) 


Pe - P5 


h - P5 

’+ 

in- 


P5 - ^4 


P6 - P5 

h|p 5 - P4 


LY5j 

(1.1.52) 


This concludes our discussion of spline vector construction. Details of 
the construction are contained in the Maintenance Document (section 4-1. 2. 3). 
We note that the "RESERVE" spline discussed there is in fact the least squares 
edge spline implemented in Pan Air. 
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Figure 1.21 - Standard half panel 



Figure 1.22 - Computation of basic far field moments 
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Neighbors for least squares fit 


Figure 1.23 Location of global source parameters used 
in the construction of the discontinuous 
source analysis spline. 
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Panel for which spline is computed 
O Location of points 


Figure 1.24 Location of source parameter points for 
the special case construction of 
discontinuous source analysis splines 
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X Edge midpoint source strengths to be 
computed by bilinear fit 
O Global source parameters 

• Global source parameters used in bilinear fits 


Figure 1.25 Computation of panel edge midpoint source 

strengths for discontinuous SD2 source splines 
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X Corner points for which spline vectors 
are to be computed 
O Global source parameters 

• Global source parameters used in bilinear fits 


Figure 1.26 Computation of panel corner spline vectors 
for continuous SD2 source splines 
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0.3 Leading and Side Edge Force 


0.3.1 Linearized Three-Dimensional Theory 

Consider a flat plate at an angle of attack, as illustrated in figure 
0.1. It is known that such a configuration experiences zero drag in subsonic 
two-dimensional potential flow, where drag is the component of force in the 
freestream or x-direction. Yet the surface normal has an x-component and the 
surface is impermeable, and thus equation (0.1.1) indicates non-zero drag. 

The resolution of this contradiction is found in the existence of a 
leading edge force resulting from an infinite leading edge velocity (see 
Ashley and Landahl, Ref. 0.1, section 5.3, or, for more detail, Hancock and 
Garner, Ref. 5.2, part II). This leading edge force exactly cancels out the 
drag computed by integrating the pressure over the configuration surface. 
Since the surface is impermeable, the momentum transfer term gives no 
contribution to the force on the surface. 

A formula for the edge force magnitude per unit distance (on a three 
dimensional wing) is given in Hancock and Garner (ref. 5.2). This formula, 
valid for all subsonic edges (there is no edge force on a supersonic edge) 
gives edge force per unit length by 


dS/dy„ = (ir/8) a [ lim (y/VT")]^ (0.3.1) 

where, S is the edge force, y„ measures distance along the edge, x„ 
measures distance perpendicular to the edge and in the plane of the surface, y 
is the doublet strength on the surface and is given by 



where Mp is the Mach number in the direction of the edge force. Thus, 

M = M„ cos A where A is the sweepback angle as shown in figure 0.2. The 
n ^ 

edge force direction is away from the edge and perpendicular to both the edge 
and the surface normal, also as shown in figure 0.2. 


0.3.2 Application in PAN AIR 

To apply equation (0.3.1) in PAN AIR, we must first consider the 
evaluation of the limit. First note that it is known that if Zp.is a 
polynomial in Xp and yp, then the exact value of y/v^ can be uniformly 
approximated by a polynomial. Given this fact, it is tempting to estimate the 
limit by evaluating y/\^ at the first and second panels from the edge and 
extrapolating these values to Xp = 0. However, because PAN AIR represents y 
itself as a quadratic polynomial in x^ and y^, the expression Mp^/ (we 

denote the doublet distribution computed by PAN AIR as yp^) will not 

converge uniformly to y/ in the neighborhood of the edge. Consequently, 
if we denote the computed and exact values of y at the center of the first 

panel by and .j^exact- "1 ,Pa'- 1, exact '■ 
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number of panels is increased. Similarly, pp^lv2 exact converge to 

1. These ratios will, however, converge to finite values, the specific values 
depending upon the solution details, especially the panel spacing and 
density. For example, different limit ratios are obtained for uniform spacing 
than for cosine spacing. Also, since a potential flow field near an edge is 
mathematically similar to two-dimensional flow about a flat plate, the 
convergence properties for three-dimensional cambered surfaces will be similar 
to those for three dimensional flat plates. Thus, we assume that the limiting 
behavior of t ® three-dimensional cambered surface will 


closely resemble the behavior of exact ^ plate as the 

number of chordwise panels is increased. 


In light of these considerations, the following method is proposed as one 
that (i) will yield reasonably accurate results without an excessive number of 
panels and (ii) will converge to exact results as the number of panels is 
increased, provided certain panel spacing rules are followed. Referring to 
figure 0.3, define Yj, Y2 by 


Y. = p/V)^ 


n,i 


i = 1,2 


(0.3.3) 


where x^ ^ denote the values of Xq at the first and second panel center 
away froili the edge (see figure 0.3 for definition). These two evaluations now 
immediately provide an extrapolated value of p/ at Xq = 0: 




(0.3.4) 


If p/ V)^ converged uniformly to its limit value we could approximate 

1 im p/ VT" = G 
x„-0 

In fact, using equation (0.3.1) we write 


dS/dy^ = (tt/8) 6^ f (0.3.5) 

where the correction factor, f (which is 0(1)) is determined from 
program/theory comparisons for a series of two-dimensional flat plate 
problems. The correction factor f depends on the number of panels and the 
panel spacing method. The total leading edge force is computed from the 
formula 


r^n,2 

S = J V (dS/dy^) dy^ (0.3.6) 

^n,l 

where the integral is evaluated by the midpoint rule using the formula (0.3.5) 
to provide values for (dS/dyp). 
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Correction factors for three types of panel spacing and various panel 
counts are given in figure 0.4. Correction factors for untabulated values of 
NPAN, the panel count, are obtained by interpolation, where f is regarded as a 
function of (1/NPAN). Correction factors have not been computed for any other 
types of panel spacing, and results obtained by PAN AIR for such types of 
panel spacing may be neither accurate nor convergent. 

The correct use and interpretation of figure 0.4 requires that one know 
precisely what is meant by "cosine spacing" and "semi-cosine spacing." Let t 
be a nondimensional coordinate defined such that the value t=0 coincides with 
the edge on which the edge force acts and t = 1 coincides with the opposite 
edge. Then the panel corner points for cosine spacing are located at the 
t-stations given by 


(cosine-spacing) t^ = (1/2)[1 - cos( (i-l)it/NPAN)] i=l, 

while for semi-cosine spacing these x stations are given by 

(semi-cosine spacing) t^- = 1 - cos[(i-l)ir/(2*NPAN)] i=l, 
0.3.3 Edge Force Verification 


.,NPAN+1 

(0.3.7) 

., NPAN+1 
(0.3.8) 


A leading edge suction distribution has been computed by PAN AIR using the 
technique described above for a 60” delta wing in incompressible flow. This 
distribution is compared in figure 0.5 to the result derived in Medan, (ref. 

0.2). Notice that the quantity (dS/dy^)/[(l-n) is plotted against n 

where n represents "spanwise fraction," the fractional transverse distance 
from the centerline to the outboard tip and a is angle of attack. This 
scaling of the leading edge suction distribution is based upon the theoretical 
results derived in (ref. 0.2). In fact, as n approaches 1, dS/dy^ itself 
tends to zero so that the differences between results shown in figure 0.5 are 
not significant when scaled and integrated to obtain overall forces. Note 
also that PAN AIR results appear to converge to the theoretical values except 
for the panels closest to the tip. In fact, the PAN AIR computed suction 
force on the tip panel will never converge no matter how densely the wing is 
panelled. This is another manifestation of nonuniform convergence arising 
from the fact that PAN AIR constrains the doublet distribution to be a 
quadratic function of the surface coordinates. Nevertheless, overall forces 
and moments do converge as panel density increases. 

Drag values predicted for a supersonic delta wing are compared in figure 
0.6 to the theoretical results presented in figure A, 14m of Jones and Cohen 
(ref. F.l). The PAN AIR results lie very close to the theoretical curve. 

Since edge suction is a significant contributor to drag, this close comparison 
provides some verification of the validity of PAN AIR'S method of computing 
edge suction forces. 
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0.5 Added Mass Coefficients 

The inertial properties of a submerged body moving in a fluid are affected 
by the motion of the fluid surrounding the body. These properties, which 
become quite significant when the mass of the displaced fluid is commensurate 
with the mass of the body, are concisely summarized by a collection of tensors 
called "added mass coefficients." In this section we will define the added 
mass coefficients, describe their computation and derive some of their 
transformation properties. 

0.5.1 Formulation of Added Mass Coefficients 

Consider an impermeable body B with boundary surface S moving through a 
fluid. We suppose that the fluid has constant density p and that the fluid 
motion is itself irrotational . Thus, a velocity potential exists and 
satisfies Laplace's equation. At the instant of observation, the motion of 
the fluid is observed from an axis system fixed in the undisturbed fluid and 
momentarily coincident with a body fixed axis system. The velocity of any 
point p on the body's surface is given by 

t + oTx (^ - fo) 

where iT specifies the body's translational velocity and If specifies its 
rotational velocity about a center of rotation Th® harmonic velocity 
potential for the fluid's instantaneous velocity is now uniquely determined by 
the impermeable surface boundary condition 

3d/an = ( iT + u X (^-po) ) • ^ on S (0.5.1) 

An equivalent form of (0.5.1) using the implied summation notation is 

3d/3n = u -j n -j to-j (^ X n ) ^ r = p ~ Pq (0.5.2) 

where we introduce the shorthand t for the vector from the center of rotation 
to a point on the body surface. Evidently, d can be written in the form 

d = u-j di ui (0.5.3) 

where di and are harmonic velocity potentials satisfying the boundary 
conditions 

3d^ /an = n. 

3 ^^/ 3 n = (f X n)^ on S (0.5.4) 

The potentials introduced in (0.5.3) are associated with motion along an axis 
in the case of di, and with rotation about an axis in the case of ♦i. 

The expression for added mass coefficients is obtained by computing T, the 
kinetic energy of the fluid induced by the motion of the body. Summarizing 
the development given in ref. 0.3, T is given by 
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T = (p/2) 

/// 

R^-B 

(V<^)^ dV 


= -(p/2) 

// 

S 

i (3(^/3n) dS 

(0.5.5) 

Substituting (0.5.3) 

into (0.5.5) yields for T 


T = (1/2) 

u. M. . 
1 ij 

“j * (1'2) “iEij Uj 


+ (1/2) 

u . S . . 
1 1J 

“j * "i ‘ij “J 

(0.5.6) 

where we define the various 

added mass coefficients 

by $he expressions 

"ij - 

S ^ 

(3^./3n) dS 

J 


S,^. = -p 

ih, 

S ^ 

(3f./3n) dS 

J 


E|j - -» 


(3(^./3n) dS 

J 


= -» 


(3^./3n) dS 

sj 

(0.5.7) 


The coefficients are called inertia ^efficients while I-jj are called 
moment of inertia coefficients; S-jj and Z/ij are called mixed'^ 
coefficients. Green's second identity implies that M and I are symmetric and 
further, that 


Thus, the full 6x6 added mass coefficient matrix given by 


“m 


E 

I 


is symmetric. 


(0.5.8) 


(0.5.9) 


0.5.2 Integration Procedures 

Because the method by which PAN AIR solves for and ensures 
internal stagnation ((#-j = = 0 interior to B), the integrals appearing 

in equation (0.5.7) can be expressed quite simply in terms of the fundamental 
singularity distributions. Letting Pi, Hi denote respectively the doublet 
distributions associated with ^i, ♦i and using the boundary conditions 
(0.5.4), we obtain 
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= - P 


"j 

dS 



^ij 

= - P 

h 

(r 

xn)j 

dS 


Sij 

= - P 

S ^ 

n . 
J 

dS 



j 

= - P 

ih, 

s ^ 

(r 

xn)j 

dS 

(0.5.10) 

Cl early these 

integrals can 

be 

evaluated by techniques 

similar to those 

described in 

section 

0.2 by 

using the 

moment matrices 

[FFmJ] and [NCPM 2 ]. 


It is important to note that the added mass coefficients given by equation 
(0.5.10) are not precisely the results computed and printed by PAN AIR. 

Rather, scaled quantities Sij, Eij and lij are printed. These 
are defined by 

"ij - 

to -Eio'Ii'Sl') 

>i0 = 

where S is a user specified reference area and L a user specified reference 
length. Notice that if one takes S = L = p = 1, one finds that = 2 


0.5.3 Orthogonal Transformation 


Added mass coefficients computed in a body axis system which has been 
subjected to an orthogonal transformation r can be defined in terms of added 
mass coefficients in a reference coordinate system. Using primes to denote 
quantities in the transformed system, the boundary conditions can be written as 


(3d-/3n') = (rn). 


(0.5.12) 


(a'l' ^ / 3 n' ) = [ r (r^ X n)]^ on S 

Equation (0.4.8) was applied to the second part of (0.5.12). For (0.5.12) to 
be satisfied, it must happen that 
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(0.5.13) 


II 


II 



Substituting (0.5.12) and (0.5.13) into a primed form of (0.5.7) will show that 
M!j - fj, (0.5.14) 

The other coefficients transform identically. 

0.5.4 Translation 

Added mass coefficients computed about an arbitrary moment reference point 
can be defined in terms of coefficients computed about the origin. The 

technique is similar to that of section 0.5.2. The boundary conditions become 

3^1 /3n' = n. 

^ 1 

3 *:/ n' = (f X n). - (f^ X n). (0.5.15) 

To satisfy these boundary conditions, the potentials in the translated system 
must take the form 

♦i - - 'ijk '■o.J \ 0-5-16) 

where denotes the usual permutation symbol, (cf. eqn. B.3.21 and ff.). 
Substitution of (0.5.15) and (0.5.16) into a primed form of (0.5.7) gives 


m; . 

ij 

= M. . 
ij 


s: . 

ij 

= S. r , M , . 

ij ikl o,k Ij 


y.'- ■ 

= y;. r , M, . 

^ij ikl o,k Ij 


I. . 
ij 

= r ,S.-E., r , y . 

IJ ilm 0,1 mj jlm o,l 



^inp ^jlm *"o,n *" 0 ,! ^pm 

(0.5.17) 


0.5.5 Symmetric Configuration 

If the surface S consists of a surface S' and its mirror image S", the 
added mass coefficients on S" can be defined in terms of those on S'. Suppose 
the plane of symmetry passes through and has a unit normal v. A point ip' 

on S' and its image point "p" on S" are related by 
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(0.5.18) 


Pi 


where 

and 


= R. . pi + X 

ij o,j 

= 6 • . - 2 V . V . 

IJ 1 J 

^o,j " ’'s.j ■ ^jk '"s,k 


ij 


The point "Xq is the image point of the origin and R is the reflection 


With primes and double primes denoting entities related to S' and S" 
respectively, the boundary conditions become 

30V/3n" = R.. nl 

34,^'/an" = - R.j (f' x + R.. (x^ x n')^ 

where t' is the vector from "p^, the center of rotation, to “p'. 


r ' = p ' - Po 

The second part of (0.5.19) required applications of (0.1.12) and (0. 
To satisfy (0.5.19), the potentials take the form 


- -"(j * «ij 


^jkl ^o,k 


Applying (0.5.19) and (0.5.20) to a form of (0.5.7) yields 


M'.' . 



S'.'. 

ij 

" ‘Jim ’'o.l "lm 

^ik ^kq *^jq 

T“- 

- -''ik^l «J1 * 

F x M" . 

ipq o,p qj 

r.'. 

IJ 

- ''ik 'kq «jq * 

F X S" 

^ikn ^o,k ^n,j 


( 0 . 


matrix. 

(0.5.19) 

5.18). 

(0.5.20) 

5.21) 
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Fp = integrated pressures 

S = leading edge force 

^ = net force 
tot 


Figure 0.1 - Effect of leading edge force 





side view 


Figure 0.2 - Leading edge force on a three-dimensional thin wing 
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Correction Factors f 
(Derived Using PAN AIR Version 1.0) 


Uniform 

Cosine 

Semi-Cosin( 

1.1579008 

1.1579008 

1.4059546 

1.3289155 

1.4039062 

1.4977210 

1.3462617 

1.4766240 

1.5680541 

1.3513135 

1.5336144 

1.5985113 

1.3544339 

1.5852960 

1.6212810 

1 .3557162 

1.6149452 

1.6323902 

1.3562257 

1.6294277 

1.6359293 

1.3564644 

1.6373295 

1.6407743 

1.3566063 

1.6431274 

1.6435113 

1.3566336 

N.A. 

1.6442279 

1.3566417 

N.A. 

1.6444954 


NPAN = Number of panel rows or number of panel columns, 
depending on which edge the edge forces are 
calculated for. 


Fi 



gure 0.4 - Edge force correction factors 
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(dS/dyf,) 






8 

Figure 0.5 - Comparison of leading edge suction predicted 
by PAN AIR with theoretical results. 


1.0 n 


0.6-4 (11/30/81) 



Figure 0.6 - Comparison of drag 
for a delta wing in 
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